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ABSTRACT 


This  paper  is  concerned  with  linear  hyperbolic  systems 
of  partial  differential  equations  for  which  certain  of  the 
associated  propagation  speeds  are  a  great  deal  larger  than  the 
other  propagation  speeds.  In  certain  cases  the  fast  modes 
allowed  by  such  a  system  are  not  present  in  the  true  physical 
s.olution.  Yet  the  fact  that  such  modes  are  allowed  means  that 
when  one  tries  to  compute  a  numerical  solution  to  an  initial¬ 
boundary  value  problem,  the  errors  generated  can  propagate 
quite  rapidly.  In  particular,  when  the  boundary  data  used 
for  the  computation  are  less  accurate  than  the  initial  data, 
the  fast  modes  can  cause  a  rapid  contamination  of  the 
calculation  in  the  interior.  To  prevent  this,  one  would 
like  to  have  boundary  conditions  which  prevent  fast  waves 
from  entering  the  region.  The  goal  of  this  paper  is  to  find 
such  conditions. 

The  situation  described  here  is  often  encountered  when 
equations  of  gas  dynamics  are  used  to  model  the  behavior  of 
the  earth's  atmosphere.  This  is  the  physical  problem  which 
motivates  this  study. 

In  order  to  find  the  desired  boundary  conditions,  we  first 
transform  the  given  system  to  an  approximate  diagonal  form  in 
such  a  way  that  each  of  the  new  dependent  variables  can  be 
identified  as  a  slow,  incoming  fast,  or  outgoing  fast  component 
of  the  solution.  We  then  find  local  boundary  conditions  which 
suppress  the  incoming  fast  part.  Pseudo-differential  operators 
are  used  throughout  the  entire  process.  The  effects  of  these 
boundary  conditions  are  analyzed  using  methods  from  the  theory 
of  propagation  of  singularities  for  linear  partial  differential 
equations . 

This  process  has  been  worked  out  in  detail  for  a  model 
problem  in  one  space  dimension  and  for  the  linearized  shallow 
water  equations,  a  system  in  two  space  dimensions.  We  have 
included  the  results  of  some  numerical  calculations  which 
demonstrate  the  effectiveness  of  the  boundary  conditions. 
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CHAPTER  1 
INTRODUCTION 

Hyperbolic  partial  differential  equations  are  characterized  by 
the  fact  that  in  a  certain  sense  they  propagate  information  at  finite 
speed.  For  first  order  hyperbolic  systems  there  may  be  several  such 
propagation  speeds,  each  corresponding  to  an  eigenvalue  of  the  princi¬ 
pal  symbol  of  the  system.  In  this  paper  we  will  consider  systems  for 
which  the  various  speeds  can  have  substantially  different  magnitudes. 
Such  systems  are  sometimes  said  to  have  "multiple  time  scales.” 

Examples  of  these  systems  arise  in  the  study  of  fluid  dynamics. 
For  such  systems  there  are  certain  propagation  modes  related  to  the 
movement  of  the  fluid,  and  there  are  certain  other  modes  which  have  a 
different  physical  interpretation.  For  the  Euler  equations  of  gas 
dynamics  these  other  modes  are  associated  with  the  movement  of  sound 
waves,  and  for  the  shallow  water  equations  they  are  related  to  the 
movement  of  gravity  waves.  If  these  waves  move  at  speeds  which  are 
considerably  greater  than  the  rate  of  flow  of  the  fluid,  then  these 
systems  have  two  time  scales. 

The  work  presented  here  is  concerned  with  a  certain  difficulty 
which  can  arise  when  one  tries  to  compute  numerical  approximat ions  to 
the  solutions  of  such  systems.  The  physical  problem  which  motivates 
this  study  is  the  use  of  hyperbolic  systems  to  model  the  behavior  of 
the  earth's  atmosphere. 
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The  specific  situation  is  illustrated  in  Figure  1.1.  This 
figure  shows  the  domain  of  definition  for  an  initial-boundary  value 
problem  for  a  hyperbolic  system  in  one  space  dimension.  In  this  case 
the  spatial  region  is  an  interval  I,  and  the  system  is  to  be  studied 
for  time  t  >  0.  The  restriction  to  one  space  dimension  is  made  solely 
for  the  purpose  of  keeping  the  picture  simple.  In  order  to  define  a 
well-posed  problem  on  this  space-time  domain,  it  is  necessary  to  specify 
values  for  the  solution  at  time  t  =  0,  and  it  is  also  necessary  to 
specify  certain  conditions  at  the  boundary  of  I  throughout  all  posi¬ 
tive  time.  In  specific  situations  the  necessary  data  are  taken  from 
physical  measurements. 


Figure  1.1 


In  certain  meteorological  problems  the  boundary  data  available 
for  a  numerical  computation  are  often  considered  to  be  substantially 
less  accurate  than  the  available  initial  data.  The  reasons  for  this 
will  be  discussed  a  little  later.  This  state  of  affairs  is  unfortunate, 
since  the  inaccuracies  in  the  boundary  data  will  generate  comparable 


inaccuracies  in  the  interior,  thereby  wasting  the  extra  accuracy 
contained  in  the  initial  data.  Our  goal  is  to  control  this  con¬ 
tamination  as  much  as  possible. 


This  sort  of  control  is  feasible  because  hyperbolic  systems 
can  propagate  information  only  at  finite  speed.  For  any  subregion 
J  of  the  given  region  I,  there  is  a  certain  period  of  time  after 
t  =  0  during  which  the  solution  in  J  cannot  be  influenced  by  the 
boundary  data.  However,  it  is  inevitable  that  the  boundary  values 
will  eventually  influence  the  solution  in  J  and  thereby  reduce  the 
accuracy  of  the  computed  solution  in  that  region.  The  question  is 
how  long  this  takes.  We  have  noted  that  the  systems  in  question  can 
allow  both  slow  and  fast  propagation  speeds.  These  are  illustrated 
by  the  characteristic  lines  appearing  in  Figure  1.1.  If  the  boundary 
data  influence  the  interior  at  the  fast  speed,  then  in  the  region  J 
the  solution  is  accurate  up  to  the  time  Tj  indicated  in  the  figure. 

If  the  boundary  data  move  in  at  the  slow  speed,  then  the  computed 
so'ution  is  accurate  up  to  time  T.,.  In  the  meteorological  problem 
these  times  can  easily  differ  by  a  factor  of  five  to  ten.  It  would 
therefore  be  worthwhile  to  prevent  this  contamination  from  taking 
place  at  the  faster  speed. 

In  order  to  do  this  we  will  try  to  find  boundary  conditions  which 
prevent  rapidly  moving  waves  from  entering  the  given  spatial  region. 

IVe  will  try  to  identify,  in  some  sense,  the  portion  of  the  solution 
which  is  entering  the  region  at  the  fast  speed,  and  we  will  then  attempt 
to  set  this  part  of  the  solution  equal  to  :ero  at  the  boundary. 
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|  For  the  meteorological  problem  this  would  accomplish  what  we 

want.  The  crucial  feature  of  this  problem  is  that  boundary  conditions 
of  this  type  are  entirely  realistic.  In  this  problem  the  fast  modes 
in  the  system  correspond  to  the  motion  of  sound  waves  or  gravity  waves. 
The  amount  of  energy  contained  in  such  waves  is  insignificant  compared 
to  the  other  forms  of  energy  in  the  atmosphere,  so  for  practical  pur¬ 
poses  the  fast  part  of  the  exact  solution  is  in  fact  equal  to  zero. 

The  only  reason  that  the  fast  modes  can  cause  any  trouble  is  that  a 
numerical  computation  can  introduce  errors  which  have  nothing  to  do 
with  the  exact  solution.  These  errors  can  therefore  be  propagated  by 
all  of  the  modes  in  the  system.  Because  the  fast  part  can  consist 
only  of  errors,  it  is  entirely  reasonable  to  try  to  suppress  this  part 
of  the  solution.  We  will  not  attempt  to  prevent  the  propagation  of 
error  at  the  slow  speed,  since  it  is  not  realistic  to  assume  that  the 
slow  part  of  the  solution  is  equal  to  zero.  We  will  instead  accept 
the  fact  that  on  any  subregion  the  computed  solution  will  eventually 
suffer  reduced  accuracy  due  to  the  effect  of  the  boundary  data. 

This  discussion  has  been  based  on  the  fact  that  there  are  certain 
modes  allowed  by  the  system  which  are  not  present  in  the  true  physical 

i 

j  solution.  These  modes  do  not  contribute  to  a  description  of  the 

I 

physical  situation,  but  they  do  cause  problems  when  we  try  to  compute 

; 

numerical  approximations  to  the  solution  of  the  system.  We  have  men¬ 
tioned  one  such  problem,  and  we  will  mention  another  a  little  later. 

It  might  seem  that  we  could  best  deal  with  these  problems  by  modifying 
the  system  of  differential  equations  so  as  to  prohibit  solutions  con¬ 
taining  rapidly  moving  waves.  This  would  certainly  eliminate  the 
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problems,  and  it  would  also  be  physically  reasonable.  In  meteorology 
this  process  is  known  as  "filtering".  However,  there  are  no  known 
filtered  systems  which  are  mathematically  well-behaved  and  which  are 
sufficiently  accurate  models  of  the  atmosphere  to  be  useful  in  meteoro¬ 
logical  calculations. 

There  is  a  partially  filtered  system,  known  as  the  "primitive 
equations",  which  is  currently  being  used  for  such  calculations.  This 
system  is  derived  from  the  Euler  equations  of  gas  dynamics,  and  it  is 
based  on  the  assumption  that  the  atmosphere  is  in  hydrostatic  balance. 
This  assumption  prevents  vertically  moving  sound  waves  from  appearing 
in  the  solution.  Unfortunately,  this  sytem  does  not  have  certain 
desirable  mathematical  properties.  The  system  is  not  hyperbolic,  and 
it  has  been  shown  by  Oliger  and  Sundstrom  [  T  ]  that  it  is  not  possible 
to  find  local  boundary  conditions  at  open  boundaries  which  lead  to 
well-posed  initial-boundary  value  problems  for  this  system.  In  current 
practice  a  diffusion  term  is  added  to  the  system  to  make  it  parabolic, 
and  values  for  all  components  are  then  prescribed  at  the  boundary. 

This  leads  to  a  well-posed  problem,  but  it  also  reduces  the  accuracy 
of  the  solution  which  is  computed. 

Because  of  the  difficulties  involved  with  finding  a  suitable 
filtered  system,  it  may  be  desirable  to  use  the  unaltered  Euler  equa¬ 
tions  of  gas  dynamics  for  meteorological  computations.  The  work  pre¬ 
sented  in  this  paper  is  concerned  with  one  of  the  difficulties  which 
can  arise  when  we  try  to  do  this. 

There  is  another  difficulty  which  can  arise  in  this  situation. 
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Because  of  the  Courant-Friedrichs-Lewv  condition,  the  fast  modes 
in  the  system  can  impose  a  severe  restriction  on  the  permissible 
time  step  for  stable  explicit  difference  approximations  to  the 
differential  equation.  In  general,  this  would  present  us  with  the 
choice  of  either  using  an  implicit  difference  method  or  an  explicit 
method  with  very  short  time  steps.  Both  choices  would  be  undesirable 
because  of  the  computational  expense  which  would  be  involved.  How¬ 
ever,  some  recent  work  by  Kreiss  [  5  ]  has  made  it  possible  to 
avoid  this  problem  by  choosing  a  suitable  set  of  initial  data.  He 
has  shown  that  certain  problems  of  instability  can  be  avoided  if  we 
smooth  our  given  data  so  that  certain  elliptic  equations  in  the 
spatial  variables  are  satisfied  at  the  initial  time.  This  makes  it 
possible  to  use  an  explicit  difference  scheme  having  a  reasonable 
time  step. 

We  still  need  to  discuss  the  reason  why  the  boundary  data  avail¬ 
able  for  certain  meteorological  computations  are  considered  to  be 
substantially  less  accurate  than  the  available  initial  data.  This 
situation  arises  in  limited  area  computations  which  are  used  to  predict 
local  atmospheric  phenomena.  Such  computations  are  made  necessary  by 
the  sice  of  the  earth's  atmosphere.  If  we  try  to  compute  the  solution 
to  a  system  of  equations  over  the  entire  atmosphere,  then  it  will  be 
necessary  to  use  an  extremely  coarse  grid  for  the  difference  equations. 
Otherwise,  the  computation  would  be  too  lengthy  for  present-day  comput¬ 
ing  machines.  In  current  practice  the  grid  spacing  for  global  atmos¬ 
pheric  computations  is  roughly  one  interval  per  two  and  a  half  degrees 


of  latitude  and  longitude.  Such  computations  can  give  useful  informa¬ 
tion  about  global  phenoman,  but  the  grid  spacing  is  too  coarse  for 
predicting  local  phenomena. 

It  is  common  practice  to  perform  additional  computations  over 
smaller  regions  with  finer  meshes  so  that  these  local  phenomena  can 
be  resolved.  For  such  a  computati  >n  the  spatial  region  is  a  cyclinder 
in  the  atmosphere  which  is  bounded  by  the  earth's  surface,  the  top  of 
the  atmosphere,  and  an  artificial  computational  boundary.  This  arti¬ 
ficial  boundary  merely  defines  the  edge  of  the  computation  and  repre¬ 
sents  nothing  physical.  The  situation  is  illustrated  schematically  in 
Figure  1.2.  The  global  computation  is  represented  by  the  coarse  grid, 
and  the  local  computation  is  represented  by  the  finer  grid  in  the  in¬ 
terval  I.  It  is  necessary  to  find  suitable  initial  data  and  boundary 

t  ;  l  I 


I  X 

Figure  1.2 

data  for  this  computation.  We  can  safely  assume  that  we  can  find 
accurate  initial  data,  since  we  would  perform  local  computations 
only  over  a  populated  region  where  there  is  a  dense  network  of  obser¬ 
vation  stations  which  are  capable  of  accurate  measurements.  The 
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problem  is  in  finding  suitable  boundary  data  when  we  are  trying  to 
predict  future  weather  patterns.  For  such  data  we  would  have  to  use 
the  results  of  our  global  prediction.  This  computation  is  made  on  a 
mesh  which  is  much  coarser  than  that  of  the  local  computation,  so  the 
results  of  this  computation  cannot  be  considered  as  accurate  as  the 
initial  data  which  are  available.  This  is  one  source  of  the  inaccuracies 
in  the  boundary  data  which  we  have  been  discussing. 

We  now  outline  the  contents  of  this  paper.  We  will  consider  initial¬ 
boundary  value  problems  for  hyperbolic  systems  having  two  time  scales. 

Our  goal  is  to  find  boundary  conditions  which  suppress  the  part  of  the 
solution  which  would  enter  the  given  spatial  domain  at  the  fast  speed. 
Although  the  systems  of  real  interest  are  quasi-linear,  we  will  consider 
only  linearized  systems.  Our  hope  is  that  a  study  of  such  systems  can 
eventually  lead  to  useful  boundary  conditions  for  the  nonlinear  problem. 

Our  basic  method  is  to  diagonalize  the  system  in  such  a  way  that 
each  of  the  new  dependent  variables  can  be  identified  as  a  slow,  in¬ 
coming  fast,  or  outgoing  fast  component  of  the  solution.  We  will  then 
attempt  to  set  the  incoming  fast  part  of  the  solution  equal  to  zero  at 
the  boundary.  Our  methods  will  rely  heavily  on  the  use  of  pseudo 
differential  operators.  In  the  Appendix  we  will  define  a  common  class 
of  such  operators,  and  we  will  state  without  proof  some  of  their  basic 
properties. 

In  Chapter  2  we  will  discuss  the  problem  for  hyperbolic  systems 
in  one  space  dimension,  and  in  Chapter  3  we  will  generalize  the  methods 
of  Chapter  2  to  problems  in  several  space  dimensions.  We  will  use  these 
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techniques  in  Chapter  4  to  derive  boundary  conditions  for  the  linearised 
shallow  water  equations.  The  results  of  some  numerical  computations  will 
be  included. 

The  work  presented  here  is  related  to  some  work  by  Engquist  and 
Majda  on  absorbing  boundary  conditions.  In  [  1  ]  they  suggested  some 
methods  for  constructing  such  conditions  both  for  scalar  wave  equations 
and  for  first  order  hyperbolic  systems.  Some  of  the  methods  which  we 
will  use  here  resemble,  in  rough  outline,  the  ideas  they  proposed  for 
hyperbolic  systems.  Their  ideas  for  scalar  wave  equations  are  developed 
in  detail  in  [  2  ]  . 
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CHAPTER  2 

THE  PROBLEM  IN  ONE  SPACE  DIMENSION 

In  this  chapter  we  consider  the  situation  for  hyperbolic  systems 
in  one  space  variable.  The  problems  of  real  interest  occur  in  more 
than  one  space  dimension,  but  certain  features  of  the  general  problem 
can  be  seen  in  this  simpler  special  case. 


2 . 1  General  Remarks 

We  will  consider  the  hyperbolic  system 


(2.1) 


for  0  <  .x  <  1,  t  >  0.  This  can  also  be  written  w  =  Aw^  +  Cw,  where 
T  2 

w  =  (u,v)  €  R  .  The  entries  in  A  and  C  are  functions  of  x  and  t. 

In  order  to  simplify  the  notation  we  have  chosen  a  system  having 
two  scalar  components.  It  can  be  seen  easily  that  the  ideas  presented 
in  this  chapter  work  equally  well  for  systems  having  several  components. 

There  is  no  loss  of  generality  in  assuming  that  A  is  diagonal. 

The  system  is  hyperbolic,  so  A  has  real  eigenvalues  and  a  complete 
set  of  eigenvectors.  If  A  is  not  diagonal,  then  a  suitable  similarity 
transformation  and  change  of  dependent  variables  can  be  made  to  bring 
the  system  to  diagonal  form. 

We  assume  |a|  <<  |b|  and  a  <  0,  b  <  0.  The  first  assumption 
guarantees  the  presence  of  propagation  speeds  having  substantially 
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different  magnitudes.  The  second  assumption  is  made  for  the  sake  of 
definiteness.  It  also  contains  the  assumption  that  det  A  ±  0,  i.e., 

that  the  boundary  is  noncharacteristic. 

The  problem  is  to  identify  the  "fast"  part  of  the  solution  of 
the  system  and  then  find  boundary  conditions  which  suppress  this  as 
much  as  possible.  To  some  degree  this  can  be  done  by  considering  the 
usual  method  of  characteristics  for  constructing  the  solution  of  the 
system.  Suppose  first  that  the  matrix  C  is  diagonal.  The  system 
(2.1)  then  uncouples  into  two  independent  equations 


au  +  c , , u 
x  11 


v 

t 


=  bv 


x 


The  first  is  an  ordinary  differential  equation  for  u  along  charac¬ 
teristic  curves  defined  bv  =  -a.  The  second  is  an  o.d.e.  for  v 

dt 

d  \ 

along  characteristic  curves  —  =  -b.  These  are  illustrated  in 

dt 

Figure  2. 1  for  the  case  ja|  <<  |bj ,  a  <  0,  b<0. 


Characteristics  for  v 


Figure  2 . 1 


12 


Initial  values  for  these  ordinary  differential  equations  are  provided 
by  initial  data  (t  =  0)  and  boundary  data  (x  =  0,  in  this  case)  for 
the  partial  differential  equations.  It  is  clear  that  data  for  v 
are  propagated  at  the  relatively  fast  speed  and  that  data  for  u 
move  at  the  slow  speed.  Boundary  conditions  which  suppress  the  fast 
part  of  the  solution  are  therefore 


u  =  given  function 
v  =  0 


x  =  0,  t  >  0  . 


Conditions  on  u  and  v  cannot  be  given  at  the  boundary  x  =  1  for 
the  example  given  here. 

The  same  boundary  conditions  work  in  the  case  where  C  is  upper 
triangular,  i.e.,  c.,^  =  0.  The  second  component  v  still  satisfies 

the  equation  v  =  bvx  *  c^v,  and  setting  v  =  0  at  the  boundary 
prevents  the  boundary  from  influencing  the  interior  at  the  fast  speed. 
In  this  case  v  appears  as  a  forcing  term  in  the  equation  for  u, 
but  this  does  not  matter  if  v  =  0. 

Trouble  can  arise  if  C  is  not  upper  triangular.  In  this  case 
u  appears  as  a  forcing  function  in  the  ordinary  differential  equations 
tor  v  along  the  characteristic  curves  —  =  -b.  Since  u  is  in 
general  nonzero,  v  will  be  nonzero  in  the  interior  even  if  it  is  set 
equal  to  zero  at  the  boundary.  The  boundary  d3ta  will  influence  the 
solution  in  the  interior  at  the  fast  speed  by  first  influencing  u, 
which  in  turn  forces  v.  The  boundary  conditions  mentioned  above  will 
of  course  have  some  desirable  effects,  but  it  would  be  better  to  have 
more  refined  boundary  conditions  which  are  more  effective  at  reducing 
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the  magnitude  of  what  propagates  in  from  the  boundary  at  the  fast 
speed. 

Such  boundary  conditions  can  be  obtained  by  transforming  the 
system  so  as  to  reduce  the  coupling  found  in  the  lower  order  term. 

This  can  be  thought  of  as  a  process  of  identifying  more  precisely  the 
quantity  which  moves  slowly  and  the  quantity  which  moves  rapidly. 
Refined  boundary  conditions  can  be  obtained  by  setting  the  new- 

fast  variable  equal  to  zero  whenever  permissible  and  then  expressing 
this  condition  in  terms  of  the  original  urknJ.  and  v. 

It  would  suffice  to  transform  C  t  ‘  i..<v  ;  riangular  form.  How¬ 
ever,  it  happens  that  when  one  uses  the  ti  .  nation  method  outlined 
in  the  next  section,  it  is  almost  as  e.;sy  to  obtain  diagonal  form  as 
it  is  to  obtain  triangular  form.  Diagonal  form  seems  a  bit  tidier,  so 
that  is  what  we  will  seek. 
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2 . 2  Outline  of  a  Method  for  Uncoupling  Systems  of  Equations 

The  uncoupling  method  used  here  is  a  method  used  by  Taylor  [10] 
to  reduce  the  coupling  in  systems  of  pseudo  differential  equations  for 
which  the  leading  symbol  is  in  block  diagonal  form.  It  is  essentially 
a  simple  perturbation  argument  which  is  disguised  by  the  language  of 
pseudo  differential  operators.  It  is  related  to  some  uncoupling 
methods  used  by  Kreiss  [  4  ]  and  O'Malley  and  Anderson  [8  ]  .  We 
will  first  outline  the  technique  using  Fourier  transforms  in  a  formal 
way.  In  later  sections  we  will  make  this  process  rigorous. 

The  svstem  is  w  =  Aw  +  Cw .  We  want  to  use  Fourier  transforms 
t  x 

to  reduce  it  to  a  system  of  ordinary  differential  equations  and  thus 
make  it  easier  to  analyze.  We  will  not  transform  in  x  because  this 
would  require  information  about  the  solution  outside  of  the  boundary. 
That  would  not  be  appropriate  in  a  discussion  of  boundary  conditions. 
Instead,  we  use  Fourier  transforms  in  t.  The  use  of  such  transforms 
will  be  justified  later  in  a  localization  argument  which  uses  proper¬ 
ties  of  pseudo  differential  operators.  These  operators  will  also  pro¬ 
vide  a  way  of  handling  equations  with  variable  coefficients.  Certain 
properties  of  these  operators  are  summarized  in  the  Appendix. 

Write  the  system  (2.1)  as 


(2.2) 


,-1  . - 

w  *  A  w  -  A  Cw  . 
x  t 


In  terms  of  components  this  is 


IS 


Introduce  formally  the  Fourier  transform  in  t.  Let  £  be  the  dual 

AAA 

variable,  and  let  u,  v,  w  be  the  transforms  of  u,  v,  w.  If  the 
coefficients  of  the  equation  are  constant,  (2.2)  becomes 

w  (x,5)  =  i£A  ^w  -  A  *Cw 
(-•3)  8  iC(A  ^  •  ]T  A  *c)w 

=  i£R(i£)w  . 

When  5  is  large  the  matrix  R(i£)  is  a  perturbation  of  the 
diagonal  matrix  A  1 .  We  will  use  a  perturbation  argument  to  reduce 
the  coupling  caused  by  the  off-diagonal  elements. 

Let  Q(i£)  =  I  +  (i S)  M,  where  M  is  a  matrix  to  be  determined. 
For  large  5,  Q  1  exists  and  has  the  expansion 

Q'1  =  I  -  rUl  *  ^(4). 

lq>  4“ 

Using  (2.5),  we  have 

QRQ'1  =  (I  +  ±  M)  (A"  1  -  i  A-  1 C)  ( I  -  iMv  (r2)) 

(2.4)' 

=  A'1  +  -X  (MA'1  -  A_1M  -  A_1n  *  r(£,'2) 

The  coupling  is  of  order  -2  if  MA  *  -  A  -  A  is  diagonal. 
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is  a  diagonal  matrix.  We  therefore  set  to  zero  the  off-diagonal  entries 
in  (2.5). 


m^b  ^  -  a  - 

3'1c12  =  ° 

(row  1 

,  column  2) 

(2.6) 

-1  .  -1 

m21a  ‘  b  m->1  ‘ 

b'lc21  =  0 

(row  2 

,  column  1) 

.4 

The  equation  (2.6)  can  be 

solved  for 

and 

m  provided  a  4  b, 

-1  i 

which  is  certainly  true  in 

this  case. 

- 

1 

-1 

a  C 12 

be  12 

< 

4 

"\ 

'  b*1  -  a'1 

a  -  b 

i 

1 

m21 

P  l21 

aC21 

-1  u'1 

a  -  b 

b  -  a 

The  entries  and  appear  only  in  the  diagonal  elements  of 

(2.6),  and  in  fact  the  terms  involving  them  cancel.  These  values  can 
therefore  be  chosen  arbitrarily.  For  convenience,  we  will  take 
m^j  =  0  and  m.,.,  =  0.  The  matrix  Q  is  then  given  by  Q  =  I  +  (i^)  ''M ,  or 
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Equation  (2.4J  then  becomes 


(2.8) 


QRQ'1  =  Q (A~ 1  -  (i4)'1A_1C)Q_1 


-1 

-a  c 


. —  1  1 
=  A  +u 
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For  later  reference  we  estimate  the  coefficient  of  the  error  term 

_  ") 

<  (£  ").  A  calculation  based  on  (2.4)  shows  that  this  coefficient  is  a 
matrix  whose  entries  are  bounded  by  the  corresponding  entries  of 


(2.9) 


constant  *  y 


Ibi*1 

i.r* 

jb'-1 

lb! -L 

Here  y  =  max{|c.^l],  and  the  constant  which  appears  first  depends 
only  on  the  number  of  terms  involved  and  the  error  made  in  approximat¬ 
ing  [ b  -  a J  1  by  [b [  1  (recall  | a J  < <  j b i ) .  In  this  case  the  con¬ 
stant  is  a  little  larger  than  4. 

We  use  (2.8)  to  reduce  the  coupling  in  equation  (2.3). 


(2.3) 


§~  (x,£)  =  USA'1  -  A*  1c)w 


(2.10) 


^  (Qw)  =  QUSA'1 


A‘1c)Q"1(Qw) 


The  entries  in  the  matrix  Q  are  independent  of  x  since  for  this 
simplified  treatment  we  have  assumed  that  the  differential  equation 


has  constant  coefficients .  Let  w  =  Qw,  and  use  the  properties  of 
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the  similarity  transformation  in  (2.8).  Equation  (2.10)  becomes 


kT.en  ^  is  large,  the  coupling  caused  by  the  lower  order  terms 

in  (2.11)  is  weaker  than  the  coupling  in  the  original  system  (2.4). 

This  means  that  for  large  £  we  can  identify  more  precisely  the  rapidly 

moving  part  of  the  solution  and  do  a  better  job  of  suppressing  it.  This 

restriction  to  high  frequencies  is  not  a  serious  one,  since  the  goal  of 

this  work  is  to  suppress  the  effect  of  numerical  errors,  which  are 

mainly  high  frequency  phenomena.  In  Section  2.6  we  will  discuss  the 

range  of  values  of  ►  for  which  the  method  is  effective. 

This  method  can  be  applied  repeatedly  to  reduce  even  further  the 

coupling  at  high  frequencies.  To  reduce  the  coupling  to  (  (£  "l,  we 

-\ 

would  multiply  equation  (2.111  by  a  matrix  of  the  form  I  +  (i^)  “M, , 
and  then  determine  M.,  in  the  same  way  th3t  we  found  the  matrix  M 
above.  In  general,  to  reduce  the  coupling  from  (  (£  ^)  to 

(\i,  n),  we  would  use  a  multiplier  of  the  form  I  +  E  nMn.  The  details 
of  this  process  involve  no  new  ideas  and  will  not  be  given  here. 

The  method  has  been  presented  for  2x2  matrices.  In  [4  ) , 

[  8  ]  ,  [  10  ]  it  is  used  for  block  matrices  having  two  square  blocks 

on  the  diagonal.  In  this  more  general  case  the  equations  correspond¬ 
ing  to  (2.6)  can  be  solved  provided  that  the  diagonal  blocks  correspond¬ 
ing  to  a  *  and  b  *  have  disjoint  spectra.  This  method  is  also 
valid  for  block  matrices  having  any  number  of  diagonal  blocks.  A 
general  form  of  this  method  will  be  discussed  in  Section  3.5. 
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We  now  find  boundary  conditions  for  this  constant  coefficient 
system  which  suppress  the  fast  part  of  the  solution  at  the  boundary. 

A  /v 

The  new  dependent  variable  is  defined  by  w  =  Ow.  By  (2.7)  this  is 


_1_  /  12  \ 

/Ul(x,C)\  /  1  iC  \  a  -  b  j  \  / u(x,5) \ 

\£(£a)  i  A*'*-55' 


v.  is  our  new  notion  of  what  constitutes  the  rapidly  moving  part  of 
the  solution.  For  large  £  it  is  a  perturbation  of  the  fast  charac¬ 
teristic  variable  v.  To  suppress  the  fast  part  of  the  solution  we 
set  v.  =  0  at  x  =  0,  i .e . , 


(2.12) 


iC  v  b  - 


+  v  =  0  at 


0  . 


To  obtain  local  boundary  conditions  we  multiply  (2.12)  by  and 

then  apply  an  inverse  Fourier  transform.  The  result  is 


3v 

3t 


aci  i 

(r~-)u  =  0 
'b-a 


at 


x 


0  . 


With  this  we  conclude  the  outline  of  the  uncoupling  method. 
There  are  several  things  left  to  do  for  problems  in  one  space  dimen¬ 
sion.  W'e  still  need  to  justify  the  use  of  Fourier  transforms  in 
time,  present  the  uncoupling  method  for  systems  having  variable  co¬ 
efficients,  and  discuss  further  the  effect  of  the  uncoupling  on  the 


behavior  of  the  solution. 


2 . 3  Transformations  in  Time 

Here  we  discuss  the  question  of  taking  Fourier  transforms  in 
time.  We  first  recall  that  it  is  a  good  idea  to  use  transforms  in 
one  variable  or  the  other  since  this  can  simplify  the  analysis  of 
the  problem.  Partial  differential  equations  can  be  reduced  to  ordinary 
differential  equations,  and  through  these  transforms  the  solutions  can 
be  expressed  as  superpositions  of  plane  waves.  This  latter  point  is 
particularly  important  for  problems  in  several  space  dimensions  since 
in  that  case  the  direction  of  propagation  can  be  as  important  as  speed. 

However,  there  are  certain  difficulties  associated  with  the  use 
of  Fourier  transforms  in  this  case.  First  of  all,  we  cannot  use 
Fourier  transforms  in  x  since  these  involve  information  about  the 
solution  outside  the  boundary.  This  is  not  appropriate  in  a  dis¬ 
cussion  of  boundary  conditions.  On  the  other  hand,  we  cannot  use 
Fourier  transforms  in  t  directly,  either.  The  reason  is  that  in  general 
the  solution  to  a  linear  hyperbolic  system  can  grow  exponentially  as 
t  -*■  +°°.  This  makes  it  impossible  to  define  a  Fourier  transform  either 
in  the  classical  integral  sense  or  in  the  sense  of  tempered  distribu¬ 
tions  . 

A  common  cure  for  this  problem  of  exponential  growth  is  the  use 
of  the  Laplace  transform.  Let  s  =  n  +  i£,  where  n  and  £  are  real 
and  n  >  0.  The  Laplace  transform  of  a  function  w  =  wft)  is  de¬ 


fined  bv 
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Lw(s)  =  !  e  St  w ( t ) dt  . 

J0 

This  is  certainly  well  defined  provided  n  is  sufficiently  large. 
However,  we  are  reluctant  to  use  this  transform  for  this  problem 
because  of  the  effect  it  has  on  the  form  of  the  differential  equation. 
Derivatives  are  transformed  according  to  the  relation 

r°° 

( Lw  1  (si  =  '  e  St  w  (t )dt 
r  J0  t 

=  s  Lw  -  w  ( 0 1  . 

The  transform  of  the  equation  w-  =  Aw  +  Cw  is  therefore 

t  x 

(2.13)  sw I x , s !  -  w(x, 0)  =  Aw  (x,s)  *  Cw  , 

x 

where  w(x,s)  is  the  Laplace  transform  in  t  for  fixed  x.  The 
trouble  with  (2.131  is  that  it  includes  initial  values  of  the  solution. 
We  would  like  to  use  our  transformed  equation  to  find  boundary  condi¬ 
tions  having  certain  properties,  but  the  presence  of  the  initial  data 
in  (2.13)  appears  to  complicate  matters. 

These  problems  with  the  Fourier  and  Laplace  transforms  can  be 
avoided  by  using  certain  properties  of  pseudo  differential  operators. 

We  were  going  to  introduce  these  operators  anyway  in  order  to  treat 
systems  with  variable  coefficients,  so  it  is  no  extra  trouble  to  use 
this  approach  to  solve  the  transformation  problem.  The  main  idea  is 
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to  localize  the  solution  in  time  in  order  to  make  the  Fourier  trans¬ 
formation  possible  and  then  show  that  this  localization  does  not  have 
a  great  effect  on  the  equation. 

We  first  recall  the  process  outlined  in  the  previous  section. 

There  we  formally  applied  a  Fourier  transformation  to  the  differential 
equation  and  then  manipulated  the  transformed  equation.  These  manipula 
tions  consisted  of  multiplying  Fourier  transforms  by  certain  functions 
of*  the  dual  variables.  In  effect,  we  were  applying  pseudo  differential 
onerators  to  both  sides  of  the  differential  equation.  The  utility  of 
these  manipulations  suffered  from  the  fact  that  the  Fourier  trans¬ 
formations  were  not  justified  and  from  the  restriction  to  systems  with 
constant  coefficients.  However,  these  problems  disappear  if  we  apply 
general  pseudo  differential  operators  directly  to  the  given  differen¬ 
tial  equation  rather  than  first  trying  to  find  a  suitable  transformed 
equation.  There  will  be  no  problem  with  variable  coefficients,  and 
the  Fourier  transformation  can  be  treated  in  the  manner  described  below 

Restrict  attention  to  a  fixed  time  interval  a  <  t  £  b,  and 
choose  e  Cq(1R)  so  that  ’JHt}  =  1  if  a  £  t  £  b.  Consider  the 
differential  equation  w^  =  Aw  +  Cw.  We  multiply  the  solution  w  by 
the  cutoff  function  t 'j  to  produce  a  function  which  has  compact  support 
in  t  and  which  therefore  has  a  Fourier  transform.  If  -;w  satis¬ 
fied  the  differential  equation,  then  in  the  case  of  constant  coeffi¬ 
cients  we  could  immediately  apply  the  transformation  to  the  differen¬ 
tial  equation.  But  this  is  not  the  case,  since  all  we  can  say  is 
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(2.14)  ~  (<I’w)  =  A  -J-  (tjiw)  +  C (ipw) 

o  i  ox 

provided  a  £  t  <_  b.  Equation  (2.14)  holds  for  t  in  this  interval 
because  =  1  there,  but  it  may  fail  to  hold  for  t  4  [a,b]. 

We  will  not  try  to  manipulate  a  transformed  equation,  but  instead 
we  will  apply  certain  pseudo  differential  operators  directly  to  the 
given  differential  equation.  In  the  next  section  we  will  construct 
operators  which  uncouple  the  equation  in  a  manner  analogous  to  that 
described  in  the  preceding  section.  The  manner  in  which  these  opera¬ 
tors  will  be  applied  requires  some  explanation. 

Write  (2.14)  as 

(2.15)  ('^'w )  ^  =  A  l(ti>w)t  -  A  *C(u'w)  . 

Denote  the  left  and  right  sides  of  (2.15)  by  L  and  R,  respectively, 
and  let  P  be  a  pseudo  differential  operator  in  t  which  we  would 
apply  to  (2. 15)  in  an  attempt  to  uncouple  the  system.  Since  L  and 
R  both  have  compact  support  in  t,  there  is  no  problem  in  forming 
P(L)  and  P(R).  The  question  is  whether  the  two  are  in  some  sense 
equal . 

We  know  that  L=R  if  a  £  t  £  b  and  L  /  R  for  certain  other 
t.  Since  pseudo  differential  operators  are  nonlocal,  we  can  conclude 
that  P(L)  and  P(R)  are  nowhere  equal  except  perhaps  at  a  few  points 
where  equality  occurs  by  accident.  But  we  can  still  say  something 
about  P(L)  -  P(R)  in  the  interval  a  £  t  £  b  where  we  know  that  L 
and  R  are  equal.  Pseudo  differential  operators  have  a  property 
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known  as  "pseudo  locality".  In  this  case  this  property  implies  that 
the  difference  P(L)  -  P(R)  on  the  interval  [a , b]  is  given  by  an 
operator  of  order  Roughly  speaking,  this  means  that  the  difference 

is  very  small  at  high  frequencies.  On  the  interval  [a,b]  we  do  not 
have  equality  of  P(L)  and  P(R),  but  instead  we  have  a  near  equality 
which  is  compatible  with  the  asymptotic  nature  of  our  method  which  was 
indicated  in  the  preceding  section. 

When  we  uncouple  the  system,  then,  we  will  first  choose  a  time 
interval  [a,b]  of  interest  and  cut  off  the  solution  outside  of  that 

interval.  This  will  make  it  possible  to  apply  pseudo  differential 
operators  to  each  side  of  the  "equation"  (2.15).  On  the  interval 
a  £  t  £  b  we  have  P(l)  =  P(R)  modulo  an  error  of  order  minus  infinity. 
This  error  term  will  be  dominated  by  other  errors  arising  in  the  un¬ 
coupling  procedure.  If  we  are  interested  in  analyzing  the  solution  on 
a  different  time  interval,  we  will  have  to  choose  a  different  cutoff 
function  '■!> .  This  will  alter  the  equations  we  obtain,  but  only  by 
modifying  the  coefficients  of  the  error  terms  in  certain  asymptotic 
formulas.  In  the  rest  of  this  paper  we  will  assume  that  the  solution 
has  been  cut  off  in  time,  and  we  will  not  bother  to  write  the  cutoff 
functions  ^  explicitly. 
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2.4  A  More  Complete  Treatment  of  the  Uncoupling  Method 


In  this  section  we  will  use  pseudo  differential  operators  to 
uncouple  the  given  system  of  partial  differential  equations.  The 
treatment  given  here  is  similar  to  that  given  in  Section  2.2,  but  it 
is  more  complete.  This  treatment  is  valid  for  systems  having  variable 
coefficients,  and  it  uses  the  method  of  Section  2.3  for  obtaining 
Fourier  transforms  in  time. 

As  before,  we  will  consider  the  equation 


(2.16) 


for  0  <  x  <  1  and  t  >  0.  We  assume  |  a |  «  !  bj  ,  and  for  the  sake  of 

definiteness  tve  will  assume  a  <  0  and  b  <  0.  The  system  can  also  be 

T  A 

written  in  the  form  w_  =  Aw  +  Cw,  where  w  =  (u,v)  e  IR“ .  The  en- 

t  x 

tries  in  A  and  C  are  functions  of  x  and  t. 

In  order  to  simplify  the  notation  we  have  chosen  a  system  having 
two  scalar  components.  The  method  presented  here  works  equally  well 
for  systems  having  several  components  and  for  systems  in  partitioned 
form. 

Write  the  system  in  the  form 

(2 .17)  w  =  A'*w  -  A_1Cw  . 

x  t 


In  terms  of  components,  this  is 
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In  Section  2.2  we  formally  applied  a  Fourier  transformation  in  t 
and  then  manipulated  the  transformed  equation.  These  manipulations 
consisted  of  multiplying  the  transformed  equation  by  matrices  of  the 
form  I  +  (i£)  Hi  and  then  determining  a  suitable  M.  This  reduced 
the  coupling  in  the  equation  to  reduce  the  coupling  from 

^  (j.  (n  D)  t0  f  (r  n}  >  we  would  use  multipliers  of  the  form 

1  *  (iCfnMn- 

In  this  section  we  will  not  try  to  manipulate  a  transformed 
equation,  for  reasons  stated  earlier.  Instead,  we  will  apply  certain 
pseudo  differential  operators  directly  to  the  given  equation  (2. 1ST. 
The  operator  which  will  reduce  the  coupling  from  order  cero  to  order 
-1  will  have  a  symbol  of  the  form  I  +  (i£)  Hi.  More  generally,  the 
operator  reducing  the  coupling  from  order  -  n+1  to  order  -n  will 
have  a  symbol  of  the  form  I  +  (i£)  nMn-  The  process  given  here  is 
similar  to  the  process  of  Section  2.2,  since  the  leading  order  term 
of  the  composition  of  two  pseudo  differential  operators  is  given  by 
the  product  of  their  symbols.  The  difference  between  this  treatment 
and  the  earlier  one  is  the  presence  of  certain  lower  order  correction 
terms  appearing  in  the  formula  for  composition  of  operators. 

Write  equation  (2. I7)  in  the  form 


(2.19) 


3w 

—  =  Gw  +  Dw  , 


where  D  =  -A_1C  and  G  is  the  operator  with  symbol  i£A  i.e.. 


Gw(x,t)  =  J  e1?t  A'^x.t)  iCw(x,C)dC 


.-1 

=  A  w^ 


Here  w(x,£)  denotes  the  Fourier  transform  of  w  with  respect  to  t 
for  fixed  x.  It  is  understood  that  w  is  cut  off  in  t  according 
to  the  remarks  of  the  preceding  section.  We  will  suppress  this  fact 
in  our  notation. 

Apply  a  pseudo  differential  operator  I+K  to  (2.19),  where  K 
is  an  operator  of  order  -1  which  is  to  be  determined.  The  symbol 
of  K  will  be  ( i^) ~ 1 M  for  some  matrix  M  depending  on  x  and  t. 
From  (2.  19)  we  obtain 


(2.20) 


~  [(I+K)w]  =  ( I*K)G(I+K) ' 1 { ( 1+K) w] 

d  X 


-l, 


(I  +  K)D(  I  +  K)  [  ( I  *  K )  w  ] 


+  K  w  . 
x 


Here  is  the  operator  whose  symbol  is  obtained  by  differentiating 

the  symbol  of  K  with  respect  to  the  parameter  x.  (I+K)  *  denotes 
a  parametrix  of  I  +  K.  It  is  not  hard  to  show  that  (I  +  K)  *  has  an 
asymptotic  expansion 


(2.21) 


(I+K)" 1  ~  I  -  K  +  K" -  KJ  +  .  .  . 


The  validity  of  this  expansion  depends  on  the  fact  that  the  order  of 
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K  is  negative. 

Let  w  =  (I+K)w.  From  (2.20)  and  (2.21)  we  obtain 

3w. 

-g~-  =  (I  +  K)G(I-K)w1  +  Dw1 

+  terms  of  order  (-1)  or  less, 
or 

3w 

(2.22)  =  Gw  +  (KG  -  GK  +  D)w 

3x  1  1 

+  (order  (-13  3 w  . 

The  operator  KG  -  GK  +  D  appearing  in  (2.221  has  order  cero.  We  want 
to  choose  K  so  that  the  leading  symbol  of  this  operator  is  diagonal, 
since  this  would  imply  that  the  coupling  in  the  system  (2.23  has  order  -1. 

Let  c  ,  c^,  denote  the  symbols  of  K  and  G,  respectively.  The 
composition  law  for  pseudo  differential  operators  implies  that  the  sym¬ 
bol  of  KG  -  GK  + D  is 


(2.23)  aKaG  -  cGaK  +  D  +  order  (-1)  . 

Let  3|.  =  (i£)  M,  where  M  is  a  matrix  depending  on  x  and  t 

which  we  shall  determine,  and  recall  that  ar  =  i£A  .  (See  (2.19).l 

b 

The  expression  in  (2.25''  then  becomes 

(2.24)  MA_l  -  A_1M  +  D  +  order  (-13  . 

This  is  the  symbol  of  the  operator  KG  -  GK  +  D  which  appears  in  (2.22). 
The  symbol  of  the  cero-order  part  of  (2.223  is  therefore  MA  *  -  A  1M  D , 
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and  the  equation  is  uncoupled  to  order  -1  if  and  only  if 

(2.25)  N1A  ^  -  A  +  D  =  a  diagonal  matrix  . 

This  is  exactly  the  conu  cion  encountered  in  Section  2.2,  equation 
(2.5).  (Recall  that  we  let  D  =  -A  in  (2.191.)  This  can  be 
solved  for  M  in  the  same  way  as  before.  From  the  earlier  work  we 
conclude  that  the  symbol  of  the  operator  K  is  given  by 


With  this  choice  of  K  the  coupling  in  equation  (2.221  has  order  -1. 

This  process  can  be  continued  indefinitely  in  order  to  reduce 

further  the  coupling  in  the  system.  To  reduce  the  coupling  to  order 

-2,  we  can  apply  an  operator  I  +  K-,  to  equation  (2.221.  where  k, 

_  -> 

has  a  symbol  of  the  form  (i^l  "M, .  The  matrix  M,  can  be  determined 
in  the  same  way  that  M  was  determined  above.  In  the  equation  for 

corresponding  to  (2.25),  the  matrix  corresponding  to  D  represents 
the  error  terms  of  order  -l  in  (2.22).  These  terms  would  have  to  be 
calculated  explicitly  when  deriving  (2.221.  The  new  dependent  variable 
for  the  system  would  have  the  form  w,  =  (I  *  K,)Wj  =  (I  +  K,)  ( I  ♦  Kj  )w. 
Further  uncoupling  can  be  carried  out  in  the  same  manner. 

We  note  that  the  symbol  in  (2.2bl  is  the  same  as  the  one  obtained 


in  Section  2.2.  This  is  not  the  case  for  symbols  which  uncouple  the 
system  further.  The  reason  for  this  is  that  these  symbols  are  influenced 
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by  error  terms  which  result  from  the  application  of  the  composition 
law  for  pseudo  differential  operators  during  prior  applications  of 
the  uncoupling  method.  These  error  terms  are  generally  nonzero  for 
systems  having  variable  coefficients ,  and  they  cannot  be  detected 
by  the  formal  treatment  of  constant  coefficient  systems  which  appeared 
in  Section  2.2. 

We  conclude  this  section  by  mentioning  a  minor  technical  diffi¬ 
culty  associated  with  a  symbol  of  the  form  (i£)  nMn .  Strictly  speak¬ 
ing,  such  a  function  cannot  be  a  symbol  of  a  pseudo  differential 
operator  because  of  the  singularity  at  ^  =  0.  But  it  can  be  modified 
in  a  neighborhood  of  E,  =  0  to  produce  a  smooth  bounded  function  of 
Such  a  change  will  affect  only  very  low  frequencies.  From  now 
on  we  will  always  assume  that  such  a  modification  h3S  been  made,  and 
we  will  ignore  this  fact  in  our  notation. 
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2 . 5  Boundary  Conditions 

We  will  now  use  the  results  of  the  preceding  section  to  find 
boundary  conditions  which  suppress  the  rapidly  moving  part  of  the 
solution.  During  the  uncoupling  process  we  adopted  a  change  of 
dependent  variables  whose  effect  was  to  weaken  the  coupling  contained 
in  the  lower  order  term,  at  least  at  high  frequencies.  The  new  depen¬ 
dent  variables  can  be  thought  of  as  more  precise  descriptions  of  the 
rapidly  moving  and  slowly  moving  parts  of  the  solution.  We  will  dinu 
our  boundary  conditions  by  attempting  to  set  the  "fast"  variable 
equal  to  :ero  at  the  boundary.  In  general,  it  is  no1-  possible  to  do 
this  exactly,  but  it  can  be  done  to  an  order  of  accuracy  which  is 
compatible  with  the  degree  of  coupling  which  remains  in  the  differen¬ 
tial  equation. 

We  will  consider  the  system  (2.22)  which  was  obtained  through 
one  application  of  the  uncoupling  method.  Systems  obtained  through 
several  applications  of  the  process  can  be  treated  in  a  similar  manner. 
Equation  (2.22)  is 

9w 

(2.27)  =  Gw  +  (diagonal  term  of  order  01  w.  +  (order  (-D)w 

I  A 

where 

(2 . 28)  w  =  ( I  +  K 1  w  . 

X 

Here  w  =  (u,v)  is  the  orginal  dependent  variable  for  the  system. 

The  symbol  of  the  operator  K  is  given  by  (2.261.  If  we  let 
w i  =  (u^/,  then  (2.26)  and  (2.28)  give 


"1 


U.29)  ux  (x,t) 


u(x,  O 


1  bC,^  .  “I 


(2.30) 


v  1  ( x ,  t ) 


ac 
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)u(x,r) 


v(x,C) 


The  components  and  v^  are  perturbations  of  the  original  com¬ 

ponents  u  and  v.  Recall  that  u  is  associated  with  the  slow 
characteristics  of  the  system  and  that  v  is  associated  with  the  fast 
characteristics.  The  new  fast  variable  is  therefore  v.,  and  our 
goal  is  to  set  this  equal  to  zero  at  the  boundary  whenever  permissible. 

This  can  be  done  easily  provided  the  system  has  coefficients 
which  are  independent  of  t.  In  this  case  the  bracketed  quantity  in 
the  integral  in  (2.30)  depends  only  on  x  and  E,  and  it  is  there¬ 
fore  the  Fourier  transform  of  v  with  respect  to  t  for  fixed  x. 

We  want  v ^  =  0  when  x  =  0.  This  can  be  accomplished  by  setting  the 
transform  equal  to  zero  at  x  =  0,  so  we  obtain 


1  1  ^  ^ 

)  U<°»S)  +  v(0,£)  =  0  for  all  >  . 

To  obtain  a  local  boundary  condition,  we  multiply  by  i 2  and  then 
invert  the  Fourier  transform.  The  result  is 


(2.32) 


3v 

at 


ac 
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b  -  a 


) u  =  0  when  x 


0  . 


This  argument  is  not  valid  if  the  coefficients  in  the  system 
depend  on  t .  In  this  case  the  bracketed  quantity  in  the  integral  in 
(2.30)  depends  on  t  as  well  as  x  and  E .  It  cannot,  therefore,  be 
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the  Fourier  transform  of  anything,  and  it  is  not  possible  to  write 
12.31'.  However,  despite  the  fact  that  the  derivation  given  above 
is  ;nvalid,  it  is  still  possible  to  use  (2.32)  as  a  boundary  condi¬ 
tion  in  the  case  of  variable  coefficients. 

Proposition  2.1.  If  u  and  v  satisfy  (2.32),  then  the  new 
"fast"  variable  satisfies 

Vj  =  (operator  of  order  -2)u  at  x  =  0  . 


Proof.  Suppose  that  (2.32)  holds,  and  apply  to  it  the  operator 
whose  symbol  is  (i£)  * .  This  gives 


i  ,,  »<  I  »♦  , 

««  '•  *  <Tr>  tb^)“  *  0 


Here  (i-)  and  "(iC)"  denote  pseudo  differential  operators 
with  symbols  (if)  1  and  if,  respectively.  The  small  circle 
denotes  composition  of  operators.  The  composition  law  for  pseudo 
differential  operators  yields 


i  ac, ,  ,,  "is  i  ^  1 ' 


+  order  (-5)  =  0  . 


The  sum  of  the  first  two  terms  is  v^,  as  can  be  seen  by  comparing 
(2.33)  and  (2.30).  The  third  term  has  order  -2.  From  this  the  result 
follows  immediately.  We  note  that  the  term  of  order  -2  is  generally 
nonzero  when  the  coefficients  of  the  system  depend  on  t. 
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The  "fast"  variable  v ^  is  therefore  small  at  high  frequencies 
when  x  =  0.  We  will  see  in  the  next  section  that  this  order  of 
accuracy  is  compatible  with  the  degree  of  coupling  remaining  in  the 
system  of  differential  equations  when  it  is  written  in  the  form  (2.27). 

We  should  note  that,  for  this  particular  case,  there  is  an  easy 
way  to  find  a  local  boundary  condition  which  sets  exactly  equal 

to  zero.  Set  (2.301  equal  to  aero,  and  write  this  as 

ac,.  t  .  _ 

v--~—  1  e 1  ^ 1  (  —  ~ )  u  d(  +  v  =  0  at  x  =  0. 

b  -  a  j  i  t. 

If  ac.^  f  0,  we  can  multiply  by  (b-al/ac^^  and  then  differentiate 
with  respect  to  t.  The  result  is 


u 


b  -  a 


ac 


)  v  ] 


21 


0  . 


The  trouble  with  this  approach  is  that  it  does  not  work  if  the  ex¬ 
pression  has  more  than  one  term  of  negative  order.  This  will  be  the 
case  if  the  system  has  three  or  more  components  or  if  we  have  applied 
the  uncoupling  method  more  than  once.  In  general,  it  is  necessary 
to  use  the  ideas  mentioned  in  the  proof  of  Proposition  2.1. 

We  will  indicate  how  this  process  works  in  the  case  where  the 
uncoupling  method  is  applied  again  to  uncouple  (2.271  to  order  -2. 

In  this  case  the  method  of  Proposition  2.1  can  be  used  to  help  generate 
a  boundary  condition,  not  just  verify  its  utility. 

After  additional  uncoupling  the  system  becomes 
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3w? 

=  (diagonal  operator  of  order  l)w7  +  (order  (-2))w 

where  w 2  =  (I  +  K-,)  (I  +  Kj)w.  Here  K2  is  a  suitably  chosen  operator 
of  order  -2,  and  =  K.  We  can  take  the  symbol  of  K7  to  be 
zero  on  the  diagonal,  as  was  the  case  for  K.  w7  can  be  written 


w2  =  (I  +  K1  +K  j  +  K2Kx)w  . 

The  operator  K7Kj  has  order  -3  and  can  be  deleted  from  the 
expression  for  the  new  dependent  variable  without  affecting  the  order 
of  the  coupling  in  the  system.  We  therefore  let 

(2. 34)  w2  =  (I  Kj  +  KJw  , 

and  the  system  becomes 
3w-, 

=  (.operator  of  order  l)w^  +  (order  (-2))w  . 
d  X  — 

-  T 

If  we  let  w.,  =  (u-,,v-,)  ,  then  v,  is  the  new  "fast"  variable. 
According  to  (2.34),  v7  can  be  written 


(2.35)  v,(x,t)» 


Cjfx.t) 

v(x,£)  *  - — -  U  + 


1- 


c,fx.t) 

(U)2 


where  c.(i£)~-  is  the  lower  left  element  of  the  svmbol  of  K. . 

I  s  I 

We  want  to  set  v?  =  0  at  x  =  0.  If  c^  and  c7  were  in- 
pendent  of  t,  then  we  could  set  to  :ero  the  bracketed  factor  in 


36 

the  integral  in  (2.35).  If  we  clear  denominators  and  invert  the 
Fourier  transform,  the  result  is 

(2.36)  — ~  +  c  ^  *  c-,u  =0  at  x  =  0  . 

)t“  1  dt 

The  method  of  Proposition  2.1  can  be  used  to  show  that  this  boundary 
condition  still  has  some  validity  when  the  coefficients  c^  and  c., 
vary  with  t.  However,  we  can  obtain  a  better  condition  for  the  case 
of  variable  coefficients  by  starting  with  a  more  general  form 

(2.3")  — V  *  c.  Ip  +  qu  =  0  at  x  =  0  , 

3t“  1  dt 

and  then  determining  q  in  a  manner  which  we  now  describe. 

3c! 

Proposition  2.2.  If  q  =  c,  -  -vp-  ,  then  the  condition  12.3") 
implies  that  the  "fast"  variable  v,  satisfies 

v,  =  (operator  of  order  -3)u  at  x  =  0  . 

If  (2.36)  is  satisfied,  i.e.  ,  q  =  c.,,  then  in  general  we  onlv  have 
v,  =  order  ( -2)  . 

Proof.  Suppose  (2.5")  holds,  and  apply  to  this  equation  the 

_  "> 

pseudo  differential  operator  whose  symbol  is  (i?)  ".  This  gives 
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«  1  »»  It  ?tt  M  1  ««  II  M 

- J  o  (iC)“  v  +  - r  o  (i£c.)  u 

(U)  (i?)“ 


It  1  tf 

♦  — o  "q”  u  =  0  . 


(U)‘ 


We  now  apply  the  composition  law  for  pseudo  differential  operators  and 
obtain 


9c, 


/K  -2  “"I  x 

+  (tt  )  u  ♦  — t  ir)  u 


It. 


U5V 


tier 


-y  u  +  (order  (-3))u  =  0 


This  simplifies  to 

"  Cj  "  "  9Cj  " 

V  +  (  JJ  )  u  +  (i£)  (q  -  2  )  u  =  (.order  ( -  S')  9  u  . 

9C1 

If  q  -  2  =  c2,  then  the  left  side  of  this  equation  is  equal  to 

vv  (See  (2.55).)  We  therefore  obtain  v7  =  order  (-5)  .  Note  that 
if  q  =  c^,  then  we  only  have  v?  =  order  (-2)  .  This  completes  the 
proof.  ■ 

From  this  it  should  be  clear  how  one  can  find  boundary  conditions 
corresponding  to  arbitrary  orders  of  uncoupling.  In  general,  when  the 
system  is  uncoupled  to  order  -n,  the  new  "fast"  variable  can  be  ex¬ 
pressed  in  terms  of  an  operator  of  order  -(n+1).  We  note  that  the 
composition  law  for  pseudo  differential  operators  can  play  an  important 
role  in  determining  these  boundary  conditions. 


i\'e  also  note  that  this  process  works  for  systems  having  more 
than  two  components.  The  dependent  variable  for  a  partially  un¬ 


coupled  system  has  the  form  w^  =  (I+K^)  •  •  (I+K^w,  where 

K.  has  order  -  j .  The  m1"^  component  of  w^  has  the  form 

(2.38)  v  +  terms  of  negative  order, 

where  v  is  the  mtfl  component  of  w.  The  process  outlined  above 
can  clearly  be  applied  to  (2.58). 

In  this  section  we  have  not  yet  discussed  boundary  conditions 
for  the  slow  part  of  the  solution.  For  the  system  (2.1)  which  we  are 
considering  here,  it  is  necessary  to  give  a  value  for  a  slow  variable 
at  x  =  0.  One  possible  condition  is 

(2.391  u  =  given  function  . 

From  (2.29)  we  can  obtain  another  condition, 

9u  i  -i 

(2.40)  —  +  (  - — -g- )  v  =  given  function. 

The  second  condition  has  little  practical  value  for  the  problem 
considered  here.  It  requires  boundary  values  for  a  derivative  of  u, 
and  in  a  numerical  computation  it  would  be  necessary  to  approximate 
this  derivative  from  measured  values  of  the  solution.  In  this  paper 
we  have  assumed  that  the  available  boundary  data  are  not  particularly 
accurate,  so  we  cannot  expect  much  accuracy  at  all  from  numerical 
differentiation.  This  implies  that  there  is  little  point  in  attempting 


Ji 
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to  use  the  boundary  condition  (2.40). 

We  therefore  propose  the  boundary  conditions 


(2.41) 


3v  /  ac21 
3t  '  b -  a 


)  u  =  0 


u  =  g,  for  x  =  0  , 


where  g  is  a  given  function  of  t.  The  first  condition  is  the  con¬ 
dition  (2.32)  which  was  obtained  through  one  application  of  the  un¬ 
coupling  method.  It  is  equivalent  to 


ft  ac0 

v (0 , t)  =  v(0,0)  +  (  7— g(t)dx  . 

J0  &  ‘  a 


The  conditions  (2.41)  thus  prescribe  values  for  the  characteristic 
variables  at  the  boundary  where  the  characteristics  enter  the  region. 
The  initial-boundary  value  problem  with  these  conditions  must  there¬ 
fore  be  well-posed.  In  Section  2.8  we  will  present  the  results  of 
some  numerical  calculations  which  compare  the  conditions  (2.41)  with 
the  simpler  conditions 

v  =  0 

u  =  given  function,  for  x  *  0  . 
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2.6  Estimates  of  the  Size  of  the  Fast  Part  of  the  Solution:  Outline 
and  Physical  Interpretation 

In  this  section  we  begin  to  examine  the  effects  of  the  boundary 
conditions  discussed  earlier.  We  will  first  estimate  the  size  of  the 
"fast"  part  of  the  solution  using  an  approach  which  has  much  the  same 
spirit  as  the  formal  treatment  of  the  uncoupling  process  given  in 
Section  2.2.  We  will  then  give  a  physical  interpretation  of  this 
result.  In  the  next  section  we  will  obtain  an  estimate  based  on  the 
more  rigorous  uncoupling  of  Section  2.4.  The  first  estimate  is  not 
rigorous  because  of  the  limitations  of  Section  2.2,  but  its  deriva¬ 
tion  is  basically  an  elementary  version  of  the  proof  of  the  second 
estimate.  We  therefore  present  the  first  in  order  to  help  explain 
and  motivate  the  other.  The  basic  method  used  here  is  essentially 
the  standard  technique  for  finding  energy  estimates  for  hyperbolic 
partial  differential  equations. 

According  to  the  discussion  in  Section  2.2,  the  system 

w^  =  Aw  +  Cw  can  be  transformed  into  a  svstem  having  the  form 
t  x  • 


(2.42) 


U,5) 


if,  A 


+  (diagonal  matrix  with  terms  of  order  zero 
or  less)  w 

n 

*  (  (5~n)w 


We  have  assumed  from  the  beginning  that  A  is  diagonal.  In  (2.42) 

w  is  given  bv 
n  " 
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«nU,C)  -  (I  *  (i^)'nMn)  *  ...  •  II  ♦  , 

where  the  matrices  are  chosen  in  the  manner  indicated  in  Section 

2.2.  Because  of  the  hypothesis  | a j  «  |b|  in  (2.1),  the  second  com¬ 
ponent  of  w^  is  the  new  "fast"  variable.  We  need  to  estimate  the 
site  of  this  component  in  solutions  which  satisfy  the  boundary  con¬ 
ditions  discussed  earlier. 

T 

Let  wn  =  (U  ,v  )  .  According  to  (2.42)  the  second  component 
satisfies  the  equation 

(2.43)  (X,S)  -  i4b_1C-n  *  Jo  (i^)-kAkvn.hn(x,-)  for  x  >  0  . 


where  hn  is  an  error  term  satisfying 


12.44) 


!hn(x- 


)  |  <  L 


—  n 


s!"n(!u(x,£) i 


I V  l  X ,  5 )  | ) 


The  functions  u  and  v  are  the  components  of  the  vector  w.  Tlie 
coefficient  b  in  (2.43)  is  taken  from  the  expression  for  the 
matrix  A  in  (2.1).  The  coefficients  \.  and  L  can  be  expressed 
in  terms  of  the  entries  in  A  and  C.  This  will  be  done  later  for 
the  case  n  =  1 . 

We  will  now  use  the  ordinary  differential  equation  (2.43)  to 

'0  i 
1  v  n  1  ‘ 

3.  If  (2.43)  and  (2.44)  hold,  then  for  x  >  0, 


estimate  the  site  of 
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|vn(x,C) i  £  !vnlO,5) |  e 


I  (R(CW  )x 


(2.45) 


+  _l_/pe2(R+a2)(x-y)  j  (yf0|2dy 
a/2\J0 


where  R(Q  =  Xq  +  (i£)  The  sum  is  taken  for  even  integers 

k.  The  constant  o  >  0  can  be  chosen  arbitrarily. 


Proof.  In  order  to  simplify  the  notation,  let  q  =  v  ,  h  =  h^, 
and  G(.C)  =  El*  ,1  (iq)  .  Equation  (2.4  5)  becomes 


(2.4b) 


q  =  if)b  q  +  G(^)q  +  h(x,^)  . 


The  subscript  denotes  differentiation. 

^  ■> 

IVe  will  find  an  inequality  for  Jqj"  and  use  this  to  estimate 
jqj".  Equation  (2.461  implies 


qxq  =  i^b_1qq  *  G(E)qq  +  h(x,£)q 


qq  =  q  (-  iCb"  1i\)  +qG(C)  q  +qh(x,q) 


Bars  denote  complex  conjugation.  The  sum  of  the  two  equations  is 


Jx  '^l2  =  +  V1 


(2. 47) 


2  Re [G(^l ]  !  q '  “  ♦  2  Re [hq] 


2R(C)  Iq|‘  *  2  Re  [  ho  ]  , 


where  R(£)  is  the  quantity  defined  in  (2.45).  The  last  equality 


follows  from  the  fact  that  the  X,  are  real.  This  in  turn  follows 

k 

from  the  derivation  of  the  partially  uncoupled  system  (2.42).  This 

derivation  is  based  upon  expansions  in  powers  of  i£,  and  all  of 

the  coefficients  in  these  expansions  are  real. 

The  second  term  on  the  right  in  (2.47)  can  be  estimated  using 
2  2 

the  inequality  2ab  <_  a  +  b  .  For  any  real  o  ^  0  we  obtain 

2  Re[hq]  £  2|  hqj  <_  2a"  \  q  I  ~  +  — t  !  hj  "  . 

2a' 

With  (2.47)  this  yields 

S  q!  2  <  2(R(0  *  a2)  |  q !  '  +  !h|:  . 

d  x  ^  - 


We  now  apply  the  Gronwall  inequality.  That  is,  we  move  the 
first  term  on  the  right  over  to  the  left  side,  multiply  by  the 
integrating  factor  exp [  -  2 (R+a~)x] ,  and  then  integrate.  The 
result  is 


! q (x)  |  _  <  e'(‘R+a  )X  ! q (O') 


+  p  e2(R.a-)(x-y)!h(>.)j:dy> 
2a“  ■'0 


To  obtain  the  final  result  (2.45)  we  recall  that  q  =  v  and  use 
the  fact  that  /a+b  <  /a  +  /b  when  a,b  £  0.  This  completes  the 
proof.  ■ 

In  the  estimate  (2.45)  there  is  a  term  which  involves  the  value 

of  v  at  x  =  0.  According  to  the  formal  treatment  of  constant- 
n 
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coefficient  problems  given  in  Section  2.2,  it  would  be  possible  to 

set  v  exact lv  equal  to  zero  at  the  boundarv  x  =  0.  In  that  case 
n  ' 

the  term  involving  v  (.0 . 2)  would  not  appear  in  (2.45").  However, 

we  saw  in  Section  2.S  that  it  is  not  quite  possible  to  satisfy 

v  =  0  when  the  coefficients  varv  with  t.  We  have  therefore  in- 
n 

eluded  the  extra  term  in  (2.451  in  order  to  suggest  the  more  general 

behavior  of  the  boundary  conditions.  According  to  Proposition  2.2 

and  the  comment  which  followed,  it  is  possible  to  obtain  v  (0,2)  = 

n 

(  (2  11  Mw.  This  is  certainly  compatible  with  the  other  term  in  (2.45), 
which  according  to  (2.44'  is  '  u.  nw) .  We  will  see  that  something 
like  this  actually  happens  with  the  more  rigorous  estimate  which  we 
will  obtain  in  the  next  section. 

The  method  used  here  to  uncouple  the  system  is  valid  asympto¬ 
tically  as  2  -  It  would  be  good  to  have  estimates  of  the  co¬ 

efficients  in  error  terms  such  as  (2.44)  in  order  to  have  a  rough 
idea  of  the  range  of  2  for  which  the  method  works.  In  the  follow¬ 
ing  proposition  we  do  this  for  (2.44)  in  the  case  n  =  1,  and  we 
also  give  the  value  of  a  relevant  parameter  appearing  in  the  estimate 
,2.  45)  . 

Proposition  2.4.  The  parameter  ,\  appearing  in  (2.45)  is 

given  bv  =  -c,,b  *.  When  n  =  1,  the  constant  L  in  v  2.44)  can 
0  22  n 

he  taken  to  be 

i  .  j 

(2.48)  L  =  constant  •  v"  jb| 

where  >  =  max  I ! c  . ! 1.  The  c..  are  the  coefficients  in  the  undifferentiated 

'  i l  1  it 
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term  in  [2.1).  The  "constant"  is  approximately  equal  to  4  in  this 
case.  (See  (2.9)  and  the  discussion  which  follows.) 


Proof.  In  (2.43)  the  parameter  \  is  defined  to  be  the 
coefficient  in  the  zero-order  term  in  the  differential  equation  for 
v  .  For  the  case  n  =  1  this  equation  is  given  by  the  second  row 
of  (2.11).  There  we  see  \n  =  -c10b  *.  This  value  does  not  change 
with  n,  since  further  uncoupling  is  obtained  through  transformations 
involving  matrices  of  the  form  I  +  (i$)  for  k  >  2.  Such  trans¬ 

formations  cannot  alter  the  term  of  order  zero. 

The  parameter  in  (2.44)  is  part  of  the  bound  on  the  error 

term  in  the  partially  uncoupled  system  (2.43''.  For  the  case  n  =  1 
this  system  is  given  by  (2.11).  The  coefficient  of  the  error  <  v 2  1 
is  bounded  by  the  matrix  in  (2.9).  From  this  the  conclusion  can  be 
read  immediately.  This  completes  the  proof.  * 


We  pause  to  interpret  this  result.  When  n  =  1  t he  "fast"  vari¬ 


able  v  satisfies 
n 


3v 


(2.49) 

where 

(2.30) 


(x,£l  =  (i£b  1  -  \iH-j  +  hjfx.^)  , 


!  h1  (x,  q)  |  £  constant  • 


*  r  I  '  ‘  I  ’  x 

uT  +  v:1 


The  first  line  is  equation  (2.45)  for  the  case  n  =  1.  The  second  is 
a  consequence  of  (2.44)  and  (2.48).  We  compare  this  to  the  situation 
in  which  we  do  not  do  any  uncoupling,  but  instead  use  the  original 
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system  w  =  Aw  +  Cw  given  in  (2.1).  In  that  case  the  "fast" 
variable  is  taken  to  be  v,  the  second  component  of  the  vector  w. 
It  satisfies 


3v 

3x 


b-\- 


- 1 

b  (c21u 


C-nVl 


The  forcing  term  in  this  equation  is  dominated  by  y|b  ^ u !  ,  since 
Y  =  max{!c^.l}.  A  comparison  with  (2. SO)  shows  that  the  uncoupling 
method  used  to  obtain  (2.49)  has  a  substantial  effect  when  '  >>y. 
This  relation  defines  what  we  will  mean  by  "large  frequencies"  in 
the  context  of  this  paper.  Ke  note  that  similar  relations  hold  more 
generally  and  that  '  and  %  both  have  the  dimensions  time 

For  large  scale  meteorological  problems  the  Coriolis  parameter 
is  usually  the  dominant  entry  in  the  coefficient  matrix  of  the  lower 
order  term  in  linearized  systems.  It  is  given  by  f  =  2.2  sin  r, 
where  2  is  the  earth's  angular  velocity  and  $  is  the  angle  of 
latitude.  That  is, 

■>-r  - 1  1  - 1 

f  =  2  •  sin  s?  hr  =  sin  hr 


The  methods  discussed  here  should  therefore  work  well  for  those  time 
frequencies  whose  order  of  magnitude  is  roughly  1  hr  *  or  greater. 

In  the  case  of  smaller  scale  problems  the  Coriolis  parameter 
may  be  dominated  by  certain  terms  which  arise  in  t he  linearization 
of  the  system.  This  will  reduce  the  maximum  wavelength  for  which 
the  uncoupling  method  is  effective.  However,  the  size  of  the  computa¬ 
tion  domain  is  also  reduced,  so  it  appea  *  that  the  method  may  still 


be  useful . 
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In  the  estimate  (.2.45)  for  v  there  are  factors  which  allow 

n 

exponential  growth  in  x  for  x  >  0.  (The  boundary  of  the  domain 
is  given  by  x  =  0.)  We  wish  to  make  a  rough  estimate  of  the  co¬ 
efficient  in  these  exponents  in  order  to  determine  the  length  scale 
on  which  this  exponential  growth  can  take  place.  One  of  these 
factors  is  exp [ (R+a“)x] ,  where 


R(£) 


n-1 


l 

k=2 


The  sum  is  taken  for  even  k,  and  the  parameter  a  can  be  chosen 
arbitrarily.  A  similar  factor  appears  in  the  integral  term  in  (2.431 
An  analysis  of  the  uncoupling  process  shows  that  the  dependence  on  k 
of  | \jJ  is  given  roughly  by  We  omit  the  details  of  this, 

but  the  main  idea  is  that  the  expansions  appearing  in  Section  2.2  are 
dominated  by  expansions  in  v '  £ '  The  behavior  of  R(£)  for  large 
%  is  therefore  governed  by  the  leading  term  A^,  which  according  to 
Proposition  2.4  is  equal  to  1 .  This  satisfies 


(2.51) 


<  JL.  = 
~  |b| 


We  recall  that  %  =  y  is  approximately  the  lowest  time  frequency  for 
which  the  uncoupling  method  can  have  an  effect.  This  corresponds  to 
a  period  of  y  *.  The  denominator  |bjy  *  in  (2.51)  is  therefore 
a  rough  approximation  to  the  length  of  the  longest  fast  wave  for  which 
the  method  applies.  This  defines  the  length  scale  on  which  the  ex¬ 
ponential  growth  can  take  place,  since  for  large  £  the  parameter  A ^ 

2 

dominates  the  coefficient  in  the  exponential  factor  exp [ (R+cO x] . 
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2 . ?  Estimates  of  the  Sice  of  the  Fast  Part  of  the  Solution: 

A  More  Complete  Treatment 

We  turn  now  to  the  problem  of  making  our  estimates  more  rigorous. 
The  estimate  (2.45)  for  the  fast  part  of  the  solution  was  obtained  by 
considering  an  ordinary  differential  equation  in  x  for  the  Fourier 
transform  in  time.  The  equation  was  obtained  through  the  uncoupling 
process  of  Section  2.2.  This  approach  gives  a  rough  idea  of  how  the 
boundary  conditions  affect  the  solution,  but  the  result  cannot  be 
considered  very  rigorous.  First  of  all,  it  obviously  cannot  work 
when  the  coefficients  in  t Ire  system  vary  with  t.  Furthermore,  this 
approach  ignores  the  problems  mentioned  earlier  regarding  Fourier  trans¬ 
forms  in  time.  A  correct  uncoupling  of  the  system  really  must  be  based 
on  the  use  of  pseudo  differential  operators,  even  if  the  system  has 
constant  coefficients,  and  a  correct  analysis  of  the  effect  of  the 
boundary  conditions  must  be  based  on  this  uncoupling.  In  this  section 
we  give  such  an  analysis. 

As  before,  we  will  obtain  estimates  which  indicate  the  behavior 
as  4  -  °°  of  the  Fourier  transform  of  the  "fast"  dependent  variable 
in  the  system.  In  the  earlier  case  we  did  this  by  estimating  v  (x,S) 
for  each  fixed  £.  However,  in  a  truly  rigorous  treatment  it  is  not 
possible  to  analyte  this  problem  one  frequency  at  a  time.  Instead, 
we  will  obtain  analogous  results  by  estimating  Sobolev  norms  of  the 
fast  part  of  the  solution.  For  any  real  number  s, 

Sobolev  space  is  defined  by 


the  norm  in  the 
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1/1 

!!  u|^  =  (  f(i  ♦  K!2)5  |u(£) |2  dO  “ 

=  !i  Asu|i  , 

L" 

where 

s/2 

(2.52)  (ASu)A  (£)  =  (1  +  U|2)  u(£)  ■ 

An  estimate  involving  a  Sobolev  norm  makes  a  statement  about 
the  behavior  of  the  Fourier  transform  as  £  -*•  From  certain  such 
estimates  we  will  be  able  to  conclude  that  the  "fast"  dependent  vari¬ 
able  is  small  at  large  frequencies.  However,  we  will  not  be  able  to 
conclude  that  the  variable  is  small  altogether,  since  the  estimates 
will  not  say  anything  about  low  frequencies.  But  this  is  not  a  severe 
loss.  The  uncoupling  process  has  an  effect  only  at  large  frequencies, 
so  it  is  only  at  such  frequencies  that  we  can  identify  the  fast  and 
slow  parts  of  the  solution.  At  low  frequencies  we  do  not  know  whether 
the  "fast"  variable  is  small,  but  on  the  other  hand  we  do  not  know 
that  it  is  really  "fast",  either. 

The  estimates  will  be  obtained  through  a  technique  which  re¬ 
sembles  the  one  used  earlier  to  obtain  (2.45).  It  is  essentially 
the  standard  technique  for  proving  energy  estimates  for  hyperbolic 
partial  differential  equations. 

According  to  the  discussion  in  Section  2.4,  the  system 

w-  =  Aw  +  Cw  can  be  transformed  into  a  system  having  the  form 
t  x 

3w 

( 2 .  53)  7t—  ( x ,  t )  =  Gw  +  2  w  +  E  w  . 

3x  n  n  n  n 
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Here  (I  is  the  operator  with  symbol 


where  a  and  b  are  defined  in  (2.1).  2^  is  a  pseudo  differen¬ 
tial  operator  in  the  time  variable.  It  has  order  zero,  and  its 
symbol  is  a  diagonal  matrix  in  which  x  may  appear  as  a  parameter. 

is  an  operator  of  order  -n  which  does  not  in  general  have  a 
diagonal  symbol  and  which  therefore  represents  the  error  in  the  un¬ 
coupling  process. 

In  (2.53)  we  should  also  include  an  error  term  which  represents 
the  effect  of  the  nrocedure  given  in  Section  2.3  for  justifying  the 
use  of  Fourier  transforms  in  time.  However,  this  term  can  be  neglected 
according  to  a  localitation  argument  which  we  will  present  a  little 
later.  iv'e  will  first  derive  the  estimates. 

T 

Proposition  2.5.  Suppose  that  (2.55)  holds,  and  let  w^ =  (u^.v^)  . 
Then  for  any  real  s  there  exist  constants  c,  and  c,  such  that  for 
x  >  0  the  "fast"  variable  v  satisfies 


,  c  x 

ii  -  .  _  i 


(2.54] 


fX  c  (x-y)  , 

+  c-,  I  e  ! | v ( y ,  • )  II"  „  dy  , 


provided  that  all  of  the  norms  are  finite.  The  norms  are  Sobolev  norms 
in  t  for  fixed  x.  The  constants  Cj  and  c,  may  depend  on  n  and 
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Here  [ A s , L ]  denotes  the  commutator  ASL  -  LAS.  We  insert  (2.57) 
into  (2.5t>)  and  obtain 

^Ikllg  =  (ASq.LAbq)  +  (LASq,ASq) 

(2- 53)  *  (ASq,[AS.L]q)  *  ([AS, L]q,ASq) 

+  CAsq,ASRw)  ♦  (ASRw,ASq)  . 

The  terms  in  (2.58)  will  be  estimated  using  the  general  fact  that 
every  pseudo  differential  operator  of  order  m  is  a  bounded  linear 
mapping  from  Ht  into  m. 

The  first  row  in  (2.58)  is  equal  to  ((L+L*)ASq,.\\i)  .  where  L* 
is  the  adjoint  of  the  operator  L.  According  to  the  Schwarz  inequalit 
this  is  bounded  by 

c.591  !|  (L+L*  ).\Sq  ||  7!|ASq|i  ,  . 

L“  I.” 

Tiie  first  factor  can  be  estimated  by  observing  that  L  +  L*  has  order 
zero,  even  though  L  and  L*  each  have  order  one.  This  is  a  con¬ 
sequence  of  the  fact  that  the  leading  symbol  of  the  adjoint  operator 

is  eaual  to  the  adjoint  of  the  leading  symbol  of  the  original  operator 

■) 

iince  L+  L*  is  therefore  a  hounded  operator  on  L",  it  follows  that 
(2.59)  can  be  bounded  by  a  constant  multiple  of  !|A  q;|"7  ,  or  !|q||“  . 

The  second  row  in  (2.581  involves  the  commutator  [A^,L]  = 

A' I.  -  LA'S.  This  operator  has  order  s,  since  the  leading  symbol  of 
the  product  of  two  operators  is  given  by  the  product  of  their  leading 
symbols.  The  commutator  is  therefore  a  bounded  mapping  from  Hh  into 


3  J 


H  (.or  L“).  It  follows  from  this  and  the  Schwarc  inequality  that 
the  second  row  of  (2.58)  can  be  bounded  by  a  constant  multiple  of 


The  third  row  in  (2.58)  is  dominated  by  a  multiple  of  j | q j U |  w  j 
since  the  operator  ASR  has  order  s-n.  From  (2.58)  we  can  there¬ 
for  conclude 


rKi 


*  K,  Uq||s!| 


s-n 


for  suitable  constants  and  Kv  This  inequality  can  be  inte¬ 

grated  in  the  same  manner  as  a  similar  inequality  which  appeared  in 
the  proof  of  Proposition  3.2.  The  result  is  (2.54).  This  completes 
the  proof.  ■ 


The  estimate  (2.54)  expresses  the  smoothness  of  the  "fast"  part 

v  in  terms  of  the  smoothness  of  its  boundary  values  v  (0,*)  and 
n  n 

j 

the  smoothness  of  the  entire  solution  w  =  (u,v)  .  According  to  the 
comments  of  Section  2.5  it  is  possible  to  choose  boundary  conditions 
for  this  one-dimensional  case  so  that  v  (0,*)  is  given  hv  an  operator 
of  order  -(n+1)  acting  on  w(0,*).  The  inequality  (2.54)  can  there¬ 
for  be  written 


c  x 

<  C_  e  1  !j  w  f  0,  -  )  I!  2 
“  J  s- n- 1 


+  c 


2  Jo 


x  c.(x-y) 

e  I,  w  (y ,  • ) 


s-n- 1 


If  w(v,«)  is  in  H  for  each  y,  then  we  can  let  s  -  n  =  r  and 

n+r 

conclude  that  v  (x,*)  is  in  H  for  each  x.  The  fast  part  of 
n  ' 

the  solution  is  therefore  n  degrees  smoother  than  the  full  solution. 
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This  argument  is  circular  as  presented,  since  the  derivation 

of  the  estimate  (2.54)  is  based  on  the  assumption  that  all  of  the 

norms  are  finite.  This  is  not  a  real  difficulty.  Suppose  that  we 

have  a  solution  w  to  (2.1)  which  satisfies  the  special  boundary 

conditions  and  which  lies  in  Hr  for  each  x.  We  need  to  show  that 

for  each  x  the  fast  part  v^  lies  in  Hr+n.  We  first  note  that  the 

r+n 

equation  (2.55),  =  Lq  +  Rw,  has  a  solution  in  H  which  is 

equal  to  v  when  x  =  0.  This  follows  from  a  little  functional 
n 

analysis  and  the  a  priori  estimate  just  obtained.  See,  for  example, 

pp.  05-65  in  Taylor  [9],  This  function  q  must  in  fact  be  equal 

to  v  for  all  x,  since  v  is  in  for  sufficients  low  s, 

n  n 

and  for  such  s  the  estimate  (2.54)  implies  uniqueness  of  solutions 
of  equation  (2.55).  We  can  conclude  that  the  fast  part  has  the  smooth¬ 
ness  properties  desired. 

One  matter  which  we  still  need  to  consider  is  the  error  term 
mentioned  earlier  which  should  have  appeared  in  (2.531.  This  term 
represents  the  effect  of  the  procedure  introduced  in  Section  2.3 
for  justifying  the  use  of  Fourier  transforms  in  time.  The  main  idea 
of  this  method  is  to  choose  a  time  interval  [a,b]  of  interest  and 
then  truncate  the  solution  outside  that  interval  by  multiplying  it  by 
a  smooth  function  with  compact  support.  This  introduces  an  error  term 
in  the  differential  equation  which,  on  the  interval  [a,b],  can  be 
represented  by  a  smoothing  operator.  Outside  [a.b]  it  cannot  in 
general  be  represented  in  such  a  manner.  This  error  term  really 
should  appear  in  (2.53),  but  it  is  possible  to  omit  this  term  if  we 
localize  the  solution  so  that  the  behavior  outside  [a.b]  becomes 


55 


irrelevant.  We  describe  how  to  do  this  now. 

The  procedure  is  based  on  a  construction  used  by  Hormander  in 
the  nroof  of  a  theorem  on  propagation  of  singularities  for  linear 
partial  differential  equations.  See  Hormander  [  3  ]  or  Nirenberg 
[  6  ]  ,  p.  44.  We  will  consider  an  equation  Pu  =  f,  where  P  is 
a  pseudo  differential  operator.  In  our  application  this  equation  is 
(2.55)  with  the  extra  error  term  added,  and  P  is  an  operator  in 
both  x  and  t.  Hormander  describes  a  method  for  localizing  the 
solution  to  a  neighborhood  of  a  given  bicharacteristic  of  P.  He 
constructs  a  pseudo  differential  operator  B  of  order  zero  so  that 
the  commutator  [ P ,  B ]  =  PB  -  BP  has  order  -  °°  and  so  that  the  s>Tnbol 
of  B  vanishes  outside  a  conical  neighborhood  of  the  bicharacterist ic . 
The  equation  Pu  =  f  implies  P(Bu)  =  BPu  +  [T,B]u,  or 

( :. oO)  P (Bu )  =  Bf  +  [P,B]u  . 

The  local  smoothness  of  u  is  given  by  the  global  smoothness  of  Bu 
(see  the  next  paragraph) ,  so  it  is  possible  to  study  the  smoothness 
of  u  along  the  bicharacteristic  by  considering  global  estimates  for 
(2.60).  This  behavior  is  not  influenced  by  the  term  [P,B]u  because 
this  term  is  automatically  smooth  in  both  x  and  t.  It  is  also  not 
influenced  by  the  values  of  f  away  from  the  bicharacteristic,  since 
these  values  are  cut  off  by  the  operator  B. 

The  fact  that  the  local  smoothness  of  u  is  determined  by  the 
global  smoothness  of  Bu  is  a  consequence  of  the  fact  that  B  is  an 
elliptic  operator  of  order  zero  in  a  neighborhood  of  the  bicharacter¬ 
istic.  To  show  this  rigorously  one  needs  to  do  a  certain  amount  of 
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work  with  cutoff  function.,.  The  necessary  arguments  are  given  in 
part  (b)  of  the  proof  of  Lemma  3  on  page  42  of  Nirenberg  [  t>  ]  . 

This  localization  process  enables  us  to  handle  the  extra  error 
term  mentioned  earlier.  Suppose  that  the  ssirbol  of  B  is  truncated 
so  that  it  is  zero  after  the  characteristic  leaves  the  time  interval 
[a,b].  Re-write  equation  (2.53)  as 

q  ( x ,  t )  =  Lq  +  Rw  +  Ew  , 

where  Ew  is  the  extra  error  term.  When  we  apply  the  operator  B 
to  this  equation,  this  term  is  replaced  by  BEw .  The  support  of  B 
in  time  is  contained  in  the  interval  [a,b] ,  and  the  singular  be¬ 
havior  of  Ew  is  confined  to  the  complement  of  [a,b].  These  facts, 
together  with  the  pseudo  local  property  of  pseudo  differential  opera¬ 
tors,  imply  that  BEw  must  be  entirely  smooth  in  t.  This  term  can 
therefore  be  treated  as  3  forcing  term  which  lies  in  any  Sobolev  class 
we  desire.  It  follows  that  the  estimate  in  Proposition  2.3  is  com- 

pletelv  valid  nrovided  that  v  is  replaced  bv  Bv  and  a  suitable 

n  f  n 

norm  of  BEw  is  inserted.  The  conclusions  about  smoothness  can  then 
be  applied  to  Bv^. 

The  method  used  here  actually  gives  more  precise  information  than 
is  implied  by  Proposition  2.5,  since  it  deals  with  propagation  along 
individual  bicharacteristics.  This  feature  will  he  useful  in  the 


study  of  problems  in  several  space  dimensions,  where  the  direction 
of  propagation  can  play  a  key  role.  In  particular,  it  will  be  more 


important  to  suppress  fast  waves  moving  in  a  direction  normal  to 
the  boundary  than  it  will  be  to  suppress  waves  moving  in  a  nearl 
tangential  direction.  The  method  given  here  allows  one  to  dis¬ 
tinguish  between  these  directions. 


In  this  section  we  present  the  results  of  some  numerical  computa¬ 
tions  involving  the  boundary  conditions  which  were  derived  earlier.  We 
consider  the  system 


C.62) 

and 


v  =  0 
u  =  g 


(x  =  0) 


(X  =  0)  . 


Here  g  is  a  given  function  of  t.  The  first  condition  in  (.2.63)  is 
the  condition  (2.32)  which  was  obtained  from  the  results  of  the  un¬ 
coupling  process. 

In  our  computation  the  system  is  approximated  by  the  leap  frog 
difference  scheme.  The  function  g  in  the  boundary  conditions  is  equal 
to  a  half  period  of  a  sine  wave  which  is  extended  by  zero.  A  for¬ 
ward  difference  is  used  to  approximate  the  derivative  in  (2.65).  The 
surfaces  pictured  in  Figures  2.5  and  2.4  are  graphs  of  ✓  u'  +  v 


as 
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a  function  of  x  and  t.  In  Figure  2.2  we  illustrate  the  con¬ 
figuration  of  these  su-Fice  plots. 


x 


Figure  2.2 

In  the  computations  we  set  the  solution  equal  to  zero  when 
t  =  0.  The  nonzero  part  of  the  solution  is  due  entirely  to  the  non¬ 
zero  boundary  data,  so  it  is  possible  to  study  the  influence  of  the 
boundary  data  by  examining  the  size  of  the  solution  in  various  parts 
of  the  (x,t)  plane.  The  solution  corresponding  to  the  simple 
boundary  condition  (2.62)  is  graphed  in  Figure  2.3,  3nd  the  solution 
corresponding  to  the  more  refined  condition  (2.b5)  is  given  in  Figure 
2.4.  it  is  clear  from  the  figures  that  the  second  condition  is  much 
more  effective  in  suppressing  the  fast  part  of  the  solution. 
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Figure  2. 


Figure  2. 


Solution  corresponding  to 
boundary  condition  (2.62). 
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CHAPTER  5 

THE  PROBLEM  IN  SEVERAL  SPACE  DIMENSIONS 

In  this  chapter  we  will  generalize  the  methods  of  the  preceding 
chapter  to  systems  in  more  than  one  space  dimension.  The  major  change 
required  lies  in  the  process  of  diagonalizing  the  leading  symbol  of 
the  differential  equation.  In  the  case  of  one  space  dimension  this 
diagonalization  causes  no  trouble  at  all,  but  in  the  case  of  several 
space  dimensions  it  can  get  rather  involved.  Otherwise,  there  is  little 
difference  between  this  case  and  the  one  discussed  earlier.  The  lower 
order  term  can  be  uncoupled  in  exactly  the  same  manner  as  before,  and 
the  fast  part  of  the  solution  can  still  be  estimated  by  localizing  the 
solution  to  a  bicharacteristic  and  then  applying  energy  estimates. 

We  will  first  discuss  the  problem  of  diagonalizing  the  leading 
symbol.  We  will  then  summarize  the  uncoupling  process  for  the  multi¬ 
dimensional  case,  and  we  will  conclude  with  a  useful  perturbation 
lemma  which  is  a  generalization  of  a  technique  that  war-  used  earlier. 

In  the  next  chapter  we  will  apply  these  metho.ls  to  the  shallow  water 
equations. 

5 . 1  Properties  of  the  Principal  Symbol 

We  will  consider  the  hyperbolic  system 

f  3 .  1 1  w  =  Aw  +  Bw  +  Cw 

t  x  y 

for  x  >  n,  y  e  IR.  Here  w(x,y,t'|  e  ]Rn,  and  A,  B,  and  C  are  real 
n  x  n  matrices  which  are  functions  of  x,  y,  and  t.  Without  loss 
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of  generality  we  will  assume  that  A  is  diagonal.  In  order  to 
simplify  the  notation  we  have  chosen  a  system  in  two  space  dimensions. 
Throughout  this  chapter  it  will  be  obvious  that  the  discussion  is 
equally  valid  for  systems  in  higher  dimensions  where  x  >  0  and 
y  e  IR^  for  k  >  2. 

There  is  no  serious  loss  of  generality  in  assuming  that  the 
spatial  domain  is  a  half-space.  If  the  given  domain  does  not  have 
this  form  but  still  has  a  smooth  boundary,  then  it  is  possible  to 
localize  the  problem  with  a  partition  of  unity  and  then  map  each 
boundary  portion  into  the  boundary  of  a  ha  If- space.  In  the  new  co¬ 
ordinates  the  problem  will  have  the  form  given  above. 

The  system  (3.1)  has  been  assumed  to  be  hyperbolic.  In  this  paper 
this  will  mean  that  for  every  real  £  and  oj  and  for  every  point 
( x , y , t ) ,  the  symbol 

(3.2)  ;A  ♦  ojB 

has  real  eigenvalues  and  a  complete  set  of  eigenvectors.  There  will 
be  no  need  to  assume  that  A  and  B  are  symmetric  or  that  the  sys¬ 
tem  is  strictly  hyperbolic. 

In  order  to  have  a  system  with  at  least  two  time  scales,  we  will 
assume  that  certain  eigenvalues  of  the  symbol  (3.2)  are  substantially 
greater  in  magnitude  than  the  others.  In  the  case  of  the  linearized 
shallow  water  equations  this  symbol  has  eigenvalues  -u  •  a  and 
-u  •  a  ±  jc|,  where  u  =  (Uj.u,,)  is  the  velocity  of  the  flow  about 
which  the  system  has  been  linearized,  and  o  is  the  vector  of  dual 
variables  (c,ta).  If  |u[  <<  c ,  then  this  sytem  has  two  time  scales. 


There  is  a  similar  set  of  eigenvalues  for  the  three-dimensional,  five- 
component  Euler  system  for  gas  dynamics.  In  this  case  the  small  eigen¬ 
value  has  multiplicity  three. 

We  now  turn  to  the  main  problem.  We  wish  to  find  boundary  condi¬ 
tions  for  (3.1)  which  prevent  rapidly  moving  waves  from  entering  the 
given  spatial  domain.  Our  plan  is  to  first  transform  the  system  to  an 
approximate  diagonal  form,  or  at  least  block  diagonal  form,  so  that 
each  of  the  new  dependent  variables  can  be  identified  as  a  slow,  incom¬ 
ing  fast,  or  outgoing  fast  portion  of  the  solution.  We  will  then 
attempt  to  set  the  incoming  fast  components  equal  to  zero. 

The  immediate  goal  is  to  diagonalize  the  leading  order  terms  in 
the  system  (3.1).  It  would  actually  suffice  to  obtain  a  block  diagonal 
form,  since  there  is  no  need  to  separate  various  incoming  fast  com¬ 
ponents  or  various  slow  components.  After  this  part  of  the  uncoupling 
has  been  accomplished,  we  can  use  the  methods  of  the  preceding  chapter 
to  reduce  the  coupling  caused  by  the  lower  order  terms. 

We  have  assumed  that  the  matrix  A  in  (3.1)  is  already  in  diagonal 
form.  This  involves  no  loss  of  generality,  since  if  A  is  not  in  that 
form  we  can  find  a  similarity  transformation  which  makes  it  diagonal  and 
then  adopt  a  suitable  change  of  dependent  variable.  Unfortunately,  it 
is  not  true  in  general  that  this  transformation  can  also  diagonalize  the 
matrix  B.  It  is  therefore  necessary  to  do  something  extra  if  we  want 
to  diagonalize  the  entire  principal  part  of  (3.1). 

In  the  case  of  constant  coefficients  it  may  be  tempting  to  use 
Fourier  transforms  in  x  and  y.  This  would  yield  the  equation 
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w(£,od,t)  =  (i£A  +  iuB) w  +  Cw  . 

The  leading  symbol  of  this  equation  can  be  diagonalized  easily  because 
it  is  a  scalar  multiple  of  the  symbol  (3.2)  discussed  earlier.  How¬ 
ever,  the  use  of  Fourier  transforms  in  x  requires  the  use  of  informa¬ 
tion  about  the  solution  away  from  the  boundary  x = 0,  and  this  is  not 
appropriate  in  a  discussion  of  boundary  conditions.  It  is  therefore 
necessary  to  take  a  different  approach. 

We  will  instead  use  Fourier  transforms  in  time  and  in  the  tangent 
variable  y.  For  the  time  being  we  will  use  these  transforms  in  a 
rather  formal  way,  and  it  will  be  understood  that  one  can  obtain  rigorous 
results  by  translating  various  arguments  into  the  language  of  pseudo 
differential  operators.  We  first  write  the  system  (3.1)  in  the  form 

(3.3)  w  =  A  -  A  *Bw  -  A  *Cw  . 

x  t  y 

It  will  be  assumed  throughout  this  discussion  that  the  matrix  A  is 

A 

invertible.  Let  w(x,uj,£)  denote  the  Fourier  transform  of  w  with 
respect  to  v  and  t  for  fixed  x.  Equation  (3.3)  implies 

(3.4)  w  (x,u>,£)  =  (i£A_  1  -  iojA^BJw  -  A_1cQ  . 

We  need  to  determine  the  values  of  oj  and  £  for  which  the  symbol 

(3.5)  £A_1  -  ojA^B 

can  be  diagonalized,  and  we  must  determine  whether  such  a  diagonal iza- 
tion  can  produce  a  transformed  system  in  which  each  component  of  the 
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dependent  variable  can  be  identified  as  slow,  incoming  fast,  or 
outgoing  fast.  The  answers  to  these  questions  are  not  immediately 
obvious,  since  we  have  chosen  a  nonstandard  set  of  variables  in  which 
to  apply  Fourier  transforms. 

In  order  to  get  started  we  must  consider  the  eigenvalues  and 
eigenvectors  of  the  symbol  (3.5).  Suppose  that  £  is  a  real  eigen¬ 
value  of  (3.5)  and  that  v  is  a  corresponding  eigenvector.  This 
means  that 

(5.6)  (CA'1  -  UjA'1B)v  =  ;v  . 

If  we  multiply  by  A  and  rearrange  the  terms,  the  result  is 
(3.  7)  (£A  +  wB)v  =  Tv  . 

The  matrix  'A  +  u>B  is  the  symbol  (5.2)  which  we  would  obtain  by  writ¬ 
ing  the  system  in  the  more  common  form  (3.1)  and  then  applying  Fourier 
transforms  in  the  usual  variables  x  and  y.  According  to  (3.6)  3nd 

(3.7) ,  this  symbol  imposes  the  same  relations  between  the  dual  vari¬ 
ables  uj,  and  t,  as  the  symbol  (3.5),  and  it  is  possible  to  find 
the  eigenvectors  of  one  symbol  by  examining  the  eigenvectors  of  the 
other.  The  difference  between  the  two  situations  is  that  in  one  case 
the  variable  £  is  treated  as  a  function  of  to  3nd  £ ,  and  in  the 
other  case  2  is  treated  as  a  function  of  £  and  to.  This  corre¬ 
spondence  between  the  two  s>"mbols  will  be  very  useful  in  studving 
(3.5).  At  this  point  in  the  discussion  we  know  a  great  deal  more 
about  (3.2)  than  we  do  about  (3.5),  and  the  correspondence  between 


the  two  will  enable  us  to  translate  information  about  one  into 
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information  about  the  other. 

We  begin  by  discussing  the  eigenvalues  of  (3.2).  In  order  to 
have  a  system  with  multiple  time  scales  we  have  assumed  that  certain 
eigenvalues  of  (3.2)  are  considerably  larger  than  the  others.  An 
example  of  such  a  set  of  eigenvalues  is  graphed  in  Figure  3.1(a). 

In  this  example  there  are  two  relatively  large  eigenvalues  and  one 
smaller  eigenvalue  for  each  £  and  w.  This  is  the  configuration  for 
the  shallow  water  equations,  and  it  is  similar  to  the  configuration 
for  the  Euler  equations  of  gas  dynamics.  In  the  latter  case  the  small 
eigenvalue  has  multiplicity  three.  Throughout  this  discussion  we  will 
assume  that  the  largest  eigenvalues  of  (3.2)  occur  in  pairs  and  have 
graphs  which  are  similar  to  the  graphs  of  the  large  eigenvalues  in 
Figure  5.1(a).  That  is,  we  will  assume  that  there  is  a  large  positive 
eigenvalue  whose  graph  is  a  narrow  cone,  though  not  necessarily  a  right 
circular  cone.  This  implies  that  there  must  also  be  a  negative  cone, 
since  if  (^,oj,^)  is  a  solution  of  (3.7)  then  so  is  (-;,  to-,  -£) . 
These  eigenvalues  generate  rapidly  moving  waves  in  all  directions.  The 
fact  that  the  graphs  are  nor  necessarily  right  circular  cones  means 
that  the  speed  can  vary  somewhat  with  the  direction  of  propagation. 

We  will  denote  by  12  the  double  cone  which  corresponds  to  the  largest 
eigenvalues,  and  we  will  denote  by  T  the  portion  of  the  (oj,0  space 
which  lies  inside  12.  These  are  labeled  in  Figure  3.1. 


(u,-')  space 


Figure  5 . 1 


He  are  now  in  a  position  to  discuss  the  eigenvalues  '  of  the 
symbol  (3.5),  1  -  uA  *B.  The  quantities  u,  and  €,  must 

satisfy  the  relation  (3.7),  which  is  the  same  as  the  relation  (3.0) 
which  was  discussed  in  the  preceding  paragraph.  We  can  therefore 
study  the  behavior  of  ;  from  graphs  like  Figure  3.1(a). 

First  of  all,  it  is  apparent  that  the  number  of  real  eigenvalues 
must  vary  with  the  position  of  (w,C) ■  If  (ai.C)  lies  in  F,  then 
there  are  two  values  of  £  which  are  associated  with  the  surface  fi  . 
One  is  positive  and  the  other  is  negative.  As  (u,£)  approaches  the 
boundary  of  F,  these  values  of  Q  approach  zero,  and  when  fu.F) 
leaves  F  the  eigenvalues  leave  the  real  axis  and  form  a  pair  of 
complex  conjugates.  The  eigenvalues  cannot  be  real,  since  for  any 
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real  C  the  point  (£,u),£)  must  lie  on  one  of  the  surfaces  in  Figure 
3.1(a).  They  are  complex  conjugates  because  they  are  eigenvalues  of 
a  real  matrix. 

It  is  safe  to  assume  that  the  other  values  of  £  do  not  behave 
in  this  manner,  at  least  in  a  neighborhood  of  F.  In  the  case  of  the 
shallow  water  equations  the  other  £  satisfies  the  equation 
£  *  u^  +  u,gj.  We  have  assumed  that  the  matrix  A  in  (5.1)  is  non¬ 
singular,  which  in  this  case  is  equivalent  to  saying  u^  /  0.  It  is 
therefore  possible  to  solve  for  £  in  terms  of  ta  and  t, ,  whether 
or  not  (cj,£)  is  in  Y.  A  similar  situation  holds  for  the  Euler  equa¬ 
tions  of  gas  dynamics.  We  will  therefore  assume  in  general  that  for 
(ui.C)  in  a  neighborhood  of  T  there  is  no  problem  in  solving  for 
the  values  of  ;  associated  with  surfaces  different  from  Q. 

We  will  now  characterize  the  behavior  of  (3.5)  when  (oo,£)  lies 
in  F.  This  is  the  only  portion  of  the  (w,£)  space  in  which  we  are 
really  interested,  since  this  is  the  only  portion  which  corresponds 
to  the  rapidly  moving  waves.  We  will  say  more  about  this  a  little 
later. 

Proposition  5.1.  If  (co , C)  r,  then  the  symbol  (3.5), 

*  -  ojA  *B,  has  real  eigenvalues  and  a  complete  set  of  real  eigen¬ 
vectors.  This  is  not  the  case  if  (oj,£)  is  not  in  F.  The  eigen¬ 
vectors  can  be  determined  from  those  of  the  symbol  (3.2),  ;A  +  wB. 

Proof.  Equations  (3.6)  and  (3.T)  show  that  the  eigenvectors  of 
(3.2)  are  also  eigenvectors  of  (3.5).  We  know  that  (3.2)  has  a  com¬ 
plete  set  of  real  eigenvectors  corresponding  to  fixed  (s>w)  and 
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various  eigenvalues  £.  IVe  want  to  show  the  same  thing  for  (3.5),  for 
fixed  (uj.C)  in  F  and  various  eigenvalues 

Suppose  that  (w,£)  is  in  F,  and  let  ^,...,5  denote  the 
eigenvalues  of  (3.5).  For  each  choose  a  basis  B.  for  the 

eigenspace  of  £^.A  +  ojB  corresponding  to  the  eigenvalue  c,.  We  are 
allowing  for  the  possibility  that  the  symbol  (3.2)  might  have  multiple 
eigenvalues.  The  elements  of  are  also  eigenvectors  of  (3.5) 

corresponding  to  the  eigenvalue  ■  We  claim  that  the  union  of  the 
B.  is  a  complete  set  of  vectors.  There  are  clearly  enough  of  these 
vectors.  The  fact  that  they  are  linearly  independent  follows  from  an 
argument  which  is  essentially  the  one  which  shows  that  eigenvectors 
corresponding  to  distinct  eigenvalues  are  linearly  independent.  Tins 
completes  the  proof.  ■ 

The  matters  discussed  in  this  section  can  be  given  a  physical 
interpretation.  Suppose  that  the  coefficients  in  (3.1)  arc  constant, 
and  let  C  =  0.  This  gives  the  system 


(3.8) 


=  Aw  +  Bw 
t  x  v 


If  we  insert  a  plane  wave  solution  v  exp(igx  +  iwy  +  iu"t)  into  (3.8), 
where  v  is  a  vector,  the  result  is  qv  =  (£A  +  u)B)v.  This  is  the 
condition  (5.7)  which  was  discussed  earlier.  The  surfaces  in  graphs 
like  Figure  3.1(a)  therefore  define  the  set  of  all  possible  frequencies 
for  plane  wave  solutions  to  (3.8).  It  is  apparent  that  the  rapidly 
moving  waves  are  associated  with  H,  which  is  why  we  are  interested 
in  the  behavior  of  (3.5)  only  for  (w,C)  in  ?. 
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In  graphs  like  Figure  3.1(a)  there  is  a  particular  wave  speed 
associated  with  each  surface  which  defines  £  as  a  function  of  w  and 
w  •  This  implies  that  it  is  possible  to  separate  fast  waves  from  slow 
waves  by  diagonalizing  the  symbol  (3.5).  It  is  also  possible  to  detect 
the  directions  in  which  the  fast  waves  are  moving.  The  portion  of  the 

surface  Q  for  which  the  product  ££  is  positive  corresponds  to 

waves  which  are  leaving  the  region  x  >  0,  and  the  portion  for  which 

<  0  corresponds  to  waves  which  are  entering  the  region.  By  pro¬ 

perly  defining  the  branches  of  £  on  the  two  sections  of  T,  we  can 
therefore  separate  the  fast  part  of  the  solution  into  incoming  and 
outgoing  components.  This  justifies  our  decision  to  seek  diagonal 
form  for  the  symbol  (5.5). 

IVe  need  to  say  a  little  more  about  the  directions  in  which  the 
various  waves  propagate.  A  plane  wave  exp(i£x  *  iwy  +  i£t)  must 
propagate  in  the  direction  i  (r,to).  If  the  point  (u>,5)  lies  on  the 
vertical  axis  in  Figure  3.1(b),  then  w  =  0,  and  the  wave  moves  in  a 
directio’  normal  to  the  boundary.  If  (u),£)  lies  near  the  edge  of 
T,  then  for  a  fast  wave  |j;|  is  small  compared  to  |oo|  ,  and  the 
wave  moves  in  a  direction  which  is  nearly  tangential.  This  observa¬ 
tion  will  be  useful  later  when  we  seek  explicit  formulas  for  bringing 
about  an  approximate  diagonal izat ion  of  the  symbol  (3.5).  The  approxi¬ 
mations  we  introduce  will  be  valid  asymptotically  as  j  ■*  0.  This 
will  lead  to  boundary  conditions  which  work  well  for  fast  waves  travel¬ 
ing  in  directions  which  have  sizeable  normal  components,  but  they  will 
not  work  well  for  waves  moving  in  directions  which  are  nearly  tangential. 


-1 

These  tangential  waves  do  not  present  any  real  problem,  since  they 
cannot  influence  the  interior  very  rapidly.  The  approximation  schemes 
are  therefore  worth  using. 
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Outline  of  the  Uncoupling  Process 


In  this  section  we  will  describe  the  uncoupling  process  for 

systems  in  more  than  one  space  dimension.  We  will  first  outline  some 

of  the  ideas  using  Fourier  transforms  in  a  formal  way,  and  we  will 

then  use  pseudo  differential  operators  to  make  the  process  more  rigorou 

We  consider  the  svstem  (3.11,  w-  =  Aw  +  Bw  +  Cw,  on  the  domain 
•  t  x  y 

x  >  0.  When  we  solve  for  w  and  apply  Fourier  transforms  in  v  and 
t,  the  result  is 

(3.9]  w  .(x,^,^J  =  (i£A  *  -  iccA  *B)w  -  A  ^Cw  . 

According  to  the  remarks  of  the  preceding  section,  the  leading  symbol 
i^A  1  -  ioiA  *B  is  diagonalizable  for  (to,£)  in  F,  and  it  is  only 
for  (j.i’i  in  F  that  we  can  have  rapidly  moving  waves.  For  the 

sake  of  neatness  we  will  use  a  cutoff  function  to  restrict  attention 

00 

to  that  set.  Let  y  be  a  C  function  of  ^  and  £  which  is  equal 
to  zero  outside  F  and  which  is  equal  to  1  on  all  of  F  except  for 
a  thin  layer  near  the  boundary.  Equation  (3.9)  can  be  written 

(3.10)  ^  =  ( i r A " 1  -  WA^Blw  -  A_1Cw  -  ( 1  -<f>)  iioA~  1  Bw  . 

The  last  term  in  (5.10)  is  an  error  term  which  is  zero  on  almost  all 
of  F.  It  is  nonzero  only  near  the  edge,  and  for  fast  waves  this 
corresponds  to  nearly  tangential  incidence.  The  error  term  is  there¬ 
fore  insignificant. 

For  any  particular  system  it  is  necessary  to  find  explicit 
formulas  for  similarity  transformations  which  bring  the  leading  symbol 
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i^A~ 1  -  wA'*B 

of  (5.10)  to  diagonal  form,  or  at  least  to  approximate  diagonal  form. 
We  note  that  it  would  actually  be  enough  to  obtain  a  block  diagonal 
form  in  which  each  block  corresponds  to  slow,  incoming  fast,  or  out¬ 
going  fast  components  of  the  solution.  This  situation  could  occur 
with  the  Euler  equations  of  gas  dynamics,  where  there  is  a  slow  mode 
of  multiplicity  three.  The  uncoupling  can  be  accomplished  either  by 
using  a  certain  perturbation  method  or  by  explicit ly  computing  *ne 
eigenvectors  of  1,3.11). 

To  use  the  perturbation  method  we  observe  that  '3.11  is  equal  to 
i  •"  times  the  matrix 

(5. 12)  A*1  -  (  f  )  ^A~ 1 B  . 

The  matrix  A  has  been  assumed  to  be  diagonal,  so  if  y  is  small, 

the  matrix  (5.12)  is  a  perturbation  of  a  diagonal  matrix.  We  can 

therefore  use  the  perturbation  argument  introduced  in  Chapter  2  to 
bring  (3.12)  closer  to  diagonal  form,  or  at  least  to  block  diagonal 
form.  We  can  applv  the  method  once  to  reduce  the  coupling  to  order 

(  v .  twice  to  reduce  it  to  order  (  ,  and  so  forth.  This  is 

■%  ^ 

one  of  the  approximation  schemes  mentioned  in  the  preceding  section 
which  work  well  in  directions  having  a  sicc-able  normal  component 
but  do  not  work  well  near  tangential  incidence.  We  will  present  a 
general  form  of  this  perturbation  method  in  the  next  section.  There 
it  will  become  apparent  that  in  the  case  of  multiple  eigenvalues  this 
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method  cannot  give  diagonal  form,  but  instead  can  give  a  satisfactory 
block  diagonal  form. 

Another  way  to  diagonalize  (3.  11)  is  to  compute  the  eigenvectors 
explicitly.  One  way  to  do  this  is  to  work  directly  with  the  matrix 
(5.11).  Another  is  to  find  the  eigenvectors  of  the  symbol  (3.2), 

'A  +  toB,  and  then  use  the  ideas  of  Proposition  3.1  to  translate  these 
vectors  into  eigenvectors  of  (3.11).  The  latter  approach  would  be 
preferable  if  (3.2)  is  easier  to  work  with  or  if  its  eigenvectors  are 
already  known. 

By  calculating  eigenvectors  we  will  be  able  to  obtain  an  exact 
diagonalization  of  (3.11)  when  (ia,£)  is  in  F.  This  may  appear  to 
be  an  advantage  over  the  perturbation  method  given  earlier.  However, 
the  expressions  for  the  eigenvectors  can  be  complicated,  and  in  order 
to  obtain  local  boundary  conditions  it  would  usually  be  necessary  to 
approximate  these  expressions  with  polynomials  or  rational  functions. 

We  would  again  use  approximations  which  are  valid  asymptotically  as 
4  0.  Although  it  does  not  give  exact  results,  the  second  approach 

allows  greater  flexibility  in  the  choice  of  approximation  methods.  The 
earlier  perturbation  approach  employs  one  fixed  method  of  approximation, 
but  here  we  have  a  choice  of  various  Taylor  approximations  or  rational 
Fade  approximations.  Engquist  and  Majda  (  [  1  )  ,  [  2]  )  found  Fade 
(approximations  particularly  useful  in  their  work  on  absorbing  boundary 
conditions  for  scalar  wave  equations. 

In  the  calculations  for  the  shallow  water  equations  which  appear 
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in  the  next  chanter,  we  will  use  the  perturbation  approach  which 
was  mentioned  first.  In  this  case  the  method  gives  sat isfactory 
results.  In  general,  however,  one  should  keep  in  mind  the  greater 
flexibility  allowed  by  the  direct  calculation  of  eigenvectors. 

We  now  give  an  outline  of  the  uncoupling  process  for  systems  in 
several  space  dimensions.  Our  intent  at  this  point  is  to  give  a 
broad  overview  of  the  method  and  avoid  details  which  could  obscure 
the  main  ideas.  We  will  go  through  the  process  in  great  detail  in 
the  next  chapter  when  we  derive  boundary  conditions  for  the  shallow 
water  equations.  These  calculations  will  be  rather  long  and  technical, 
so  it  will  be  worthwhile  to  first  see  a  relatively  short  outline  of 
the  process. 

We  first  solve  for  w  in  (3.1)  to  obtain  the  form  (3.5). 

w  =  A~\v  -  A"  lBw  -  A ' 1 C w  . 
x  t  y 

In  order  to  simplify  the  notation  we  will  change  the  meaning  of  A, 

B,  and  C  and  write  the  system  as 

( 3 . 13)  w  =  Aw  +  Bw ^  *  Cw  . 

A,  B,  and  C  will  have  this  meaning  throughout  the  remainder  of 
this  paper.  The  matrix  A  is  diagonal. 

In  order  to  prepare  for  the  uncoupling,  ve  will  express  (3.15) 


in  terms  of  certain  pseudo  differential  operators.  As  in  Chapter  2 
the  solution  w  will  be  truncated  in  t  so  that  these  operators 
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can  be  applied  to  it,  but  this  fact  will  be  suppressed  from  the  nota¬ 
tion.  Denote  by  H  the  operator  with  symbol 

(3.  14 J  <jh  =  i^A  +  iw^B  , 

and  let  E  be  the  operator  with  symbol  iwCi-vOB.  The  system  (3.13) 
then  becomes 

(3.15)  w  =  Hw  +  Cw  +  E,w  . 

x  1 

As  before,  is  a  smooth  function  which  is  equal  to  1  on  almost  all 
of  r  and  is  equal  to  zero  outside  T .  The  operator  represents 

an  error  which  is  zero  when  <£  =  1.  In  the  case  of  variable  coeffi¬ 
cients  the  cone  H,  and  therefore  T,  can  vary  with  x,  y,  and  t, 
so  £  is  in  general  a  function  of  x,  y,  t,  oj,  and  E,. 

We  have  mentioned  that  when  we  try  to  uncouple  the  leading  symbol 
it  will  be  useful  to  use  perturbation  arguments  involving  the  quantity 
y  .  The  operators  which  transform  the  system  must  therefore  involve 
this  quantity.  A  potential  problem  with  this  is  that  cannot  be  the 
symbol  of  a  pseudo  differential  operator  because  of  the  singularity  in 
the  direction  £  =  0.  However,  we  have  avoided  this  difficulty  by  our 
use  of  the  function  IVe  will  find  that  the  ratio  j  can  appear 

only  in  the  form  ^  <p,  and  this  is  nonzero  only  in  a  conical  neighbor¬ 
hood  of  the  axis  u  =  0.  The  singularity  is  thereby  eliminated. 

When  we  uncouple  the  system  (5.15)  the  first  task  is  to  take  care 
of  the  leading  order  operator  H.  Let  q  be  a  matrix  such  that 
qo^q  ^  is  anproximately  diagonal  or  approximately  block  diagonal, 


and  let  Q  be  the  pseudo  differential  operator  whose  symbol  is  q. 

The  operator  Q  must  have  order  zero,  since  its  symbol  q  is  homo¬ 
geneous  of  order  zero  in  its  dependence  on  ta  and  £.  If  we  apply 
Q  to  (3.15),  the  result  is 

(3.16)  (Qw)x  =  (QHQ'^Qw  +  (QCQ'1  +  QxQ_1)Qw  +  QEjW  . 

Here  Q  *  denotes  a  parametrix,  or  approximate  inverse,  of  Q.  This 
operator  is  defined  by  the  property  that  QQ  *  -  I  is  a  smoothing 
operator.  It  is  not  hard  to  show  that  such  an  operator  exists  and 
to  obtain  an  asymptotic  expansion  for  its  symbol.  An  outline  of 
the  argument  is  given  in  Section  4.1.  The  leading  order  term  in 
the  expansion  is  q  *,  the  inverse  of  the  symbol  of  Q. 

We  need  to  examine  the  operator  QHQ  *.  According  to  the  composi 
tion  law  for  pseudo  differential  operators,  its  leading  symbol  is  the 
product  of  the  leading  symbols  of  Q,  H,  and  Q  This  is  the  matri 
qa(1q  *  which  is  known  to  be  approximately  uncoupled.  There  are  also 
various  lower  order  terms  in  the  expansion  of  QHQ  These  are  due 

partly  to  the  effect  of  the  composition  law  and  partly  to  the  lower 
order  terms  in  the  expansion  for  Q  The  composition  law  is  stated 
in  the  Appendix. 

Let  G  be  the  pseudo  differential  operator  whose  symbol  is  the 
diagonal  part,  or  block  diagonal  part,  of  qo^q  \  and  let  R  be 
the  operator  whose  symbol  is  the  rest  of  qa^q  *.  G  and  R  both 
have  order  one.  Because  of  the  approximate  uncoupling  in  q-t^q  1 , 
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the  effect  of  R  is  small  except  near  tangential  incidence.  If  we 
let  wQ  =  Qw,  (3.16)  becomes 

3w 

(3-17)  TF=  Gw0  +  Zw0  +  Rwo  +  E2wo- 

Here  Z  is  the  pseudo  differential  operator  associated  with  the  zero- 
order  terms  appearing  explicitly  in  (3.16)  and  with  the  terms  of  order 
zero  or  less  which  arise  in  the  expansion  of  QHO  1 .  E is  equal  to 
QE^Q  and  its  symbol  is  equal  to  a  smoothing  term  when  =  1 .  The 

system  (3.17)  is  uncoupled  near  normal  incidence,  up  to  terms  of  order 
zero. 

The  coupling  in  the  lower  order  term  can  be  reduced  by  using  the 
same  technique  that  was  used  in  Chapter  2  for  systems  in  one  space 
dimension.  That  is,  we  can  apply  to  (3.17)  an  operator  of  the  form 
I  +  K,  where  1  is  the  identity  operator  and  K  is  an  operator  of 
order  -1  which  is  to  be  determined.  When  we  apply  the  operator 
I+K,  the  result  is 

(3.18)  ~  [ (I+K)wq]  =  (I+K)G(I+K)'1w1 

>  (I  +  K) (Z+R+E^) (I+K) _1w^  +  Kxwq  , 

where  =  (I+K)w^  =  'I+K)Qw.  The  parametrix  (I+K)  1  has  the 

a 

asymptotic  expansion  I  -  K +  K“  -  ...  .  This  follows  easily  from  the 
fact  that  the  order  of  K  is  negative.  The  system  (3.18)  can  there¬ 


fore  be  written 
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3w 

=  GWj  +  (KG  -  GK  +  Z)w 

(3.19)  +  (terms  of  order  -1  or  less) 

+  (K+K) (R+E2) CI+K) _ 1w1  . 

The  zero-order  coupling  in  (3.19)  is  caused  by  the  operator 
KG  -  GK  +  Z.  Its  leading  symbol  is 

(3.2°)  aKaG  -  aGaK  +  “0  ’ 


where  and  are  the  symbols  of  K  and  G,  and  z^  is  the 

leading  symbol  of  Z.  In  order  to  eliminate  the  coupling  of  order 
zero,  we  will  need  to  determine  so  that  the  symbol  (3.20)  is 

diagonal  (or  block  diagonal).  In  Chapter  2  we  did  this  calculation 
for  a  special  case,  and  in  the  next  section  we  will  give  a  more 
general  treatment  as  part  of  a  general  perturbation  lemma.  There  we 
will  find  that  it  is  possible  to  find  a  suitable  provided  that 

the  diagonal  blocks  of  have  disjoint  spectra.  This  condition  can 
be  satisfied  here  since  we  are  trying  to  separate  slow,  incoming  fast, 
and  outgoing  fast  components  of  the  solution. 

The  technique  given  here  can  be  used  to  uncouple  the  system 

further.  To  reduce  the  coupling  from  order  -n+  1  to  order  -n,  we 

would  use  an  operator  of  the  form  I  +  K  ,  where  K  has  order  -n. 
f  n  n 

After  m  uncouplings  the  dependent  variable  would  be 


=  (I  +  K  ) 
m 


(I+K  )Qw  , 
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where  is  the  operator  K  discussed  above. 

Boundary  conditions  for  the  system  can  be  generated  and  then 
analyzed  using  ideas  which  are  similar  to  those  used  in  the  case 
of  one  space  dimension.  To  do  the  analysis  we  would  localize  the 
solution  to  a  neighborhood  of  a  bicharacteristic  and  then  find  energy 
estimates  involving  Sobolev  norms.  These  estimates  would  give  informa¬ 
tion  about  the  behavior  of  the  solution  at  high  frequencies.  The  fact 
that  the  solution  can  be  localized  to  a  bicharacteristic  means  that 
we  can  study  the  effects  of  the  boundary  conditions  for  various  angles 
of  incidence  to  the  boundary.  This  gives  meaning  to  the  use  of  the 
approximations  about  normal  incidence  which  were  mentioned  earlier. 
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3. 5  A  Perturbation  Lemma 

In  this  section  we  present  a  method  for  reducing  the  coupling 
found  in  matrices  which  are  perturbations  of  block  diagonal  matrices. 
This  method  can  be  used  to  partially  uncouple  the  leading  symbol  in 
the  system  (3.13),  and  it  is  essentially  the  method  which  has  already 
been  used  to  reduce  the  coupling  caused  by  lower  order  terms.  We  pre¬ 
sent  it  as  a  separate  lemma  for  the  sake  of  clarity  and  generality. 
Various  versions  of  this  method  have  been  used  in  [4],  [8],  and 

[10]  for  2x2  block  matrices. 

Proposition  3.2.  Let  A  and  B  be  square  matrices  of  equal 
dimension.  Suppose  that-  A  is  block  diagonal,  and  let  Aj,...,A 
denote  the  blocks  on  the  diagonal.  If  no  two  of  the  A^  have  any 
eigenvalues  in  common,  then  for  small  s  the  sum  A  +  eB  can  be 
uncounted  to  order  e".  More  precisely,  there  exists  a  matrix  M 
such  that  for  z  sufficiently  small, 

(I  +  eM)  (A+eB)  (I+eM)"1  =  A  +  e  •  (block  diagonal  matrix)  +  (  (O  . 

A  method  for  constructing  M  will  be  given  in  the  proof. 

Proof.  For  small  z  the  inverse  ( I*eM)  *  exists  and  is  equal 

*)  'f 

to  I  -  rM  +  e“M“  -  ...  .  We  can  therefore  write 

CI  +  eM)  (A*eB)  (I+gM)'1  =  (I+EM)  (A*eB)  (I-eM  +  (  (s2) ) 

=  A  +  G(MA-AM+B)  ♦  r{ Z~)  . 
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Our  goal  is  to  choose  M  so  that 
(5.21)  MA  -  AM  +  B 

is  block  diagonal.  For  the  sake  of  notation  we  will  partition  M  and 

B  into  block  structures  which  match  the  block  structure  of  A.  M. . 

il 

and  B^.  will  denote  the  blocks  in  the  (i,j)  position.  They  are 

not  necessarily  square,  since  we  are  not  assuming  that  and  Aj 

have  the  same  dimensions.  The  (i,j)  block  in  (3.21)  can  then  be 

written  as  M.  .A.  -  A.M.  +  B.  ..  For  i  ji  i  .  we  want  this  to  be  equal 
i.i  J  i  lj  ij 

to  cero.  We  are  therefore  faced  with  the  problem  of  solving  the 
equation 


(3.22) 


M. .A.  -  A.M.  .  =  -  B.  . 
lj  1  i  11  il 


for  M...  Once  we  have  done  this,  the  proof  is  complete.  There  are 
no  conditions  imposed  on  the  diagonal  blocks  ,  so  these  may  be 
chosen  arbitrarily. 

If  A.  and  A.  are  both  1  x  1  matrices,  i.e.,  scalars,  then 
il 

we  obviously  need  to  have  A^  t  Aj  in  order  to  be  able  to  solve  (3.22) 
for  arbitrary  B._..  In  the  general  case  the  system  (3.22)  is  solvable 
if  and  only  if  A^  and  Aj  have  disjoint  spectra.  Proofs  of  this 
fact  can  be  found  in  several  different  references.  We  give  one  here 
for  the  sake  of  completeness. 

In  order  to  simplify  the  notation  we  will  write  (3.22)  in  the  form 
XS  -  TX  =  Y,  where  S,  T,  and  Y  are  given  and  X  is  to  be  deter¬ 
mined.  We  are  assuming  that  S  and  T  are  square  matrices  which  do 
not  have  any  eigenvalues  in  common.  There  is  no  need  to  assume  that 
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they  have  the  same  dimension.  We  will  denote  the  columns  of  X  and 

Y  by  x.  and  y.  ,  and  we  will  denote  the  entries  of  S  bv  s... 

.11  ij 

We  can  assume  that  S  is  upper  triangular,  since  otherwise  we 
can  use  a  similarity  transformation  to  reduce  the  problem  to  that  case. 
We  will  solve  for  the  columns  of  X,  starting  from  the  left.  We  first 
have  -  Tx^  =  y^.  The  matrix  s^I  -  T  is  nonsingular  since 

s^  is  an  eigenvalue  of  S  and  therefore  not  an  eigenvalue  of  T. 

The  column  is  therefore  determined  uniquely.  We  next  have 

-  Tx0  =  v0  -  s  ,Xj.  This  system  has  a  unique  solution  x, 
since  s10  is  not  an  eigenvalue  of  T.  We  can  continue  in  this 
manner  to  solve  for  all  of  X.  We  note  that  the  condition  on  the 
eigenvalues  of  S  and  T  is  necessary  as  well  as  sufficient,  since 
it  is  equivalent  to  the  statement  that  s^I  -  T  is  nonsingular  for 
all  i.  This  completes  the  proof  of  the  lemma,  and  therefore  the 
proof  of  the  main  proposition. 


CHAPTER  4 


AN  EXAMPLE  IN  TWO  SPACE  DIMENSIONS 

In  this  chapter  we  will  use  the  methods  of  Chapter  3  to  derive 
boundary  conditions  for  the  linearized  shallow  water  equations.  The 
calculations  will  follow  the  outline  given  in  Section  3.2.  When  we 
uncouple  the  leading  symbol  of  the  system  we  will  use  the  perturba¬ 
tion  method  given  in  Section  3.3.  This  process  and  the  one  used  to 
reduce  the  lower  order  coupling  will  each  be  applied  one  time.  Cer¬ 
tain  portions  of  the  calculations  are  specific  to  the  shallow  water 
equations,  but  other  portions  are  more  generally  applicable.  For 
much  of  the  chapter  the  spatial  domain  we  consider  will  be  the  half¬ 
space  x  >  0,  but  later  we  will  discuss  the  effect  of  rotation  of 
coordinates  on  the  form  of  the  boundary  conditions.  In  the  last  sec¬ 
tion  we  will  present  the  results  of  some  numerical  tests  of  these 
conditions . 

4 . 1  Uncoupling  the  System 

The  linearized  shallow  water  equations  can  be  written  in  the  form 
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In  this  notation,  (-a,-S)  is  the  velocity  of  the  flow  about 
which  the  system  has  been  linearized,  and  c  is  a  speed  associated 
with  the  propagation  of  gravity  waves.  IVe  will  assume  ja, , j  3 i  <<  c 
and  a  i  0.  The  dependent  variables  are  given  by  u  =  cu',  v  =  cv', 
and  p  =  (J)',  where  u'  and  v'  are  the  perturbations  in  the  com¬ 
ponents  of  velocity  and  is  the  perturbation  in  the  geopotential. 

The  coefficients  •  in  the  undifferentiated  term  are  due  partly  to 
Coriolis  effects  and  partly  to  the  process  of  linearization.  For  the 
time  being  we  will  consider  the  system  on  the  domain  x  >  0. 

The  first  step  is  to  diagonalize  the  coefficient  matrix  of  the 
normal  derivative  —  in  (4.1).  To  do  this  we  use  the  matrix 

0  3  1  v 

/2  0  0  J  . 

0  1  -1  / 


M 


The  columns  of  this  matrix  are  normalized  eigenvectors  of  the  co¬ 
efficient  matrix  in  question.  When  we  multiply  (4.1)  on  the  left  by 
the  inverse  of  (4.2)  and  make  the  appropriate  change  of  dependent 
variable,  the  result  is  the  system 


(4.3) 


a-c 


+  Dw 


a+c 


where 
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D  and  E  can  be  written  explicitly,  but  we  will  wait  until  we  solve 

for  w  . 
x 

IVe  can  use  (4.3)  to  make  a  preliminary  identification  of  the 
slow,  incoming  fast,  and  outgoing  fast  parts  of  the  solution.  Since 
|a|  <<  c,  we  can  say  roughly  that  the  first  component  v  is  a  slow 
component  and  that  u  +  p  and  u  -  p  are  incoming  and  outgoing  fast 
components,  respectively.  In  order  to  suppress  the  incoming  fast 
part  of  the  solution  we  could  therefore  require  u  +  p  =  0  at  the 
boundary  x  =  0.  If  a  <  0,  then  we  would  also  prescribe  a  value 
for  v  in  order  to  have  a  well-posed  problem.  The  trouble  with  this 
approach  is  that  it  ignores  the  effect  of  the  terms  Dw^  and  Ew. 

Our  identification  of  the  various  parts  of  the  solution  is  therefore 
not  very  accurate.  The  purpose  of  the  uncoupling  process  is  to  pro¬ 
duce  a  more  accurate  identification  and  thereby  enable  us  to  find 
boundary  conditions  which  are  more  effective  at  suppressing  the  in¬ 
coming  fast  part. 

We  need  to  solve  for  w  in  (4.5).  When  we  do  this  the  result 

x 

is 

w  =  Aw  +  Bw  +  Cw  , 
x  t  y 

where 


j^ie. 
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qcuq  ^  is  closer  to  diagonal  form  than  c,,,  at  least  for  small 
H  H 

To  do  this  we  will  use  the  ideas  of  sections  3.2  and  3.3  to 
find  a  matrix  M  such  that 


The  off-diagonal  elements  in  (4.101  are  determined  uniquely  by  the 
condition  (4.9),  but  the  diagonal  elements  may  be  chosen  arbitrarily. 
For  convenience  we  have  set  these  equal  to  cero. 
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We  now  define  the  symbol  q  by 

(4.12)  q  =  I  ♦  , 

where  M  is  given  in  (4.10),  and  we  let  Q  be  the  pseudo  differen 
“ial  operator  with  symbol  q.  When  the  operator  Q  is  applied  to 
the  system  (4.7),  the  result  is 

(4.15)  (Qw)x  =  (QHQ_1)Qw  ♦  (QCQ'1  *  QxQ_1)Qw  ♦  QEjW  . 

Here  denotes  the  operator  with  symbol  q^. 

_  1 

We  saw  in  Section  5.2  that  the  leading  symbol  of  OHQ  1  is 
qc^q  * .  According  to  (4.8)  and  (4.12)  this  is  given  by 

qCMq  1  =  (I  +  y<£M)(i£A  +  iu.'^B)(I  *  y  1  . 

If  we  factor  out  i£  and  identify  with  z  in  (4.  IP,  we  can 

A  2 

conclude  that  qc^q  1  is  equal  to  +  {%-•#'),  where 
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We  will  let  G  denote  the  operator  whose  symbol  is  0^.  The  system 
(.4.13)  can  then  be  written  in  the  form 

3w  n  2  7 

(4.15)  =  Gwq  +  ZWq  +  C(  +  ^(w(W))w  , 

where  w =  Qw. 

The  operator  Z  is  associated  with  the  zero-order  terms  appearing 

explicitly  in  (4.13)  and  with  the  terms  of  order  zero  or  less  which 

arise  in  the  expansion  of  the  symbol  of  QHQ  *.  In  a  moment  we  will 

u2  .  ' 

discuss  this  further.  The  term  f'(  —  ^“)w  denotes  the  effect  of  an 

W  -> 

operator  whose  symbol  is  dominated  by  and  which  is  a  result 

of  the  error  in  uncoupling  the  leading  order  part  of  the  system  (4.7). 
The  term  (  (aj(l-^))w  represents  the  effect  of  the  operator  Ej  which 
appears  in  (4.13).  Its  symbol  is  equal  to  zero  on  almost  all  of  the 
set  7  ,  which  is  the  only  part  of  the  (o),Q  space  in  which  fast 
waves  can  be  found.  In  F  it  is  nonzero  only  near  the  edge,  which 
for  fast  waves  corresponds  to  nearly  tangential  incidence.  This  term 
is  therefore  of  no  consequence. 

The  system  is  now  partially  uncoupled  near  normal  incidence,  since 
the  symbol  of  G  is  diagonal.  We  next  need  to  reduce  the  coupling 
caused  by  the  zero-order  operator  Z,  which  is  given  by 

Z  =  QCQ  1  +  Q^Q’1  *  (terms  order  zero  or  less 
arising  from  the  expansion  of  QHQ  *)  . 


(4.16) 


The  coupling  can  be  reduced  by  the  method  presented  in  Chapters 
2  and  3,  but  we  will  first  have  to  identify  the  leading  symbol  of 


We  first  consider  the  first  two  terms  in  (4.16).  The  leading 
symbol  of  QCQ  1  is 

qCq'1  =  (I  ♦  ♦  |  ^M) 

=  C  +  ('{  |  <fi)  , 


and  the  leading  symbol  of  Q^Q 


is 


q  q" 1  -  (7*MJ(I  ♦  y  <£M  1  ” 1 


= 


The  expression  for  q  is  taken  from  (4.121 .  We  will  regard  terms 
of  order  <f  as  error  terms  since  ^  <p  is  no  larger  than  the 

(i)**  ^ 

term  -r-  ^  which  has  already  appeared  in  (4.15).  The  first  two 
terms  in  (4.16)  are  then  given  by 


(4.17)  QCQ1  +  QxQ_1  =  C  *  order  (-1)  *  (  ( |  v?) 

In  order  to  consider  the  expansion  of  QHQ  1  we  must  first 
find  the  symbol  of  the  parametrix  Q  1 .  We  start  by  assuming  that 
the  symbol  has  an  asymptotic  expansion  of  the  form  rQ  ■  , 

where  r^  is  homogeneous  of  order  -k  in  u  and  £.  The  ciioice 
of  orders  will  turn  out  to  be  appropriate,  since  Q  has  order  zero 


L 
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We  will  then  solve  for  the  r.  one  bv  one. 

k 

Let  be  the  operator  with  symbol  r^.  Accord;ng  to  the 

composition  law  for  pseudo  differential  operators,  the  symbol  of 


1  3q  3r0  1 
i  3?  3t  i 


3oi 


3r0 

3y 


+  order  (~2) 


The  composition  law  is  stated  in  the  Appendix.  If  we  choose 
r^  =  q  then  the  leading  term  is  I.  We  note  that  r^  really 
is  of  order  zero  and  that  the  error  term  really  is  of  order  -2. 
We  now  have 


(4.18) 


QRq  =  I  +  order  (-1)  . 


Now  choose  r^  so  that 

1  3rn  3rn 

(4.19)  qrt  =  -  T(qcir*  %  ^  ■ 

The  leading  symbol  of  QR^  therefore  cancels  the  leading  symbol  of 
the  error  in  (4.18).  This  gives 

Q (Rg  +  Rj )  =  I  +  order  (-2)  . 


This  process  can  be  continued  indefinitely  to  find  any  r^.  At  each 
step  we  would  need  to  know  the  leading  symbol  of  the  error  in  the 
equation 
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This  can  be  calculated  from  the  expansions  of  the  symbols  of  QR. 
for  j  <  k. 

In  the  present  problem  we  really  only  need  the  terms  r0  and 
r j .  Q  has  order  zero  and  H  has  order  one,  so  the  terms  r,,r,,... 
must  contribute  terms  of  negative  order  in  the  expansion  of  QHQ 
These  lower  order  terms  are  of  no  interest  to  us  here. 

The  leading  symbol  Tq  of  Q  *  is  given  by 


(4 .20)  rn  =  q'1  =  (I  *  £  tpM)  =  I  *  <  (  t  *?)  • 

u  s 

We  will  see  later  that  in  this  case  it  is  not  necessary  to  keep  any 
more  terms  in  this  expansion  of  r^.  The  second  term  r^  in  Q  ' 
is  defined  in  (4. 19)  to  be 


r!  =  iq 


( qr 


?ro 

3t 


3r0 

?v 


)  ■ 


Equations  (4.12)  and  (4.20)  can  be  used  to  obtain 

r,  =  i(I  +  ^Vll'V  --^*M)  ♦  (  i  *M)  <  (  f  *) } 

L  s  r-  s  ^  s 

(4.21) 

=  f'i  ^V) 

5“ 

In  (4.21)  we  have  omitted  derivatives  of  v  since  these  are  nonzero 


only  near  the  edge  of  T  and  cannot  be  of  any  consequence.  From  the 
above  work  we  can  conclude  that  the  symbol  of  Q  *  has  the  expansion 
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The  symbol  of  QHQ  *  is  given  by 


1  3q  3 

QHQ  1  HQ  1  1  at  HQ  1 


*  f  IS  w°  -i  *  ord"(-1’  • 

HQ 


We  now  use  (4.22)  and  the  fact  that  q  is  equal  to  1  + 


-1 


a  ,  =  q [crHq  «■  (  (?  *)  +  order  (-1)  ] 
QHQ  1  H  s 


(4.23) 


+  1  (-  “  *M)  0  , 

1  3t  HQ  1 


+  f  (  j  y?M)  4r.  c  ,  *  order  ( - 1 ) 
1  5  3y  HQ'1 


A  short  calculation  shows  that 


_L 

3t 


(i£A.  +  iw^B  )  ( I  -  £*M) 

t  Z  ^ 


+  (i£A  +  iw*B)(-  T  >pM. ) 

S  t 


2 

♦  *  order  (-1). 


The  same  relation  holds  when  we  replace  t  with  v.  The  second  term 

in  (4.25)  is  therefore  ('{  ^  >f>)  *  order  (-1),  and  the  third  is 

>£MAV  +■  f  ( £  >p)  +  order  (-1).  Equation  (4.23)  can  then  be  simplified  to 
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(4.24)  O  .  =  qo-.q’1  +  tfMA  +  (’{  f  v)  +  order  (-1)  . 

QHQ'  Y  t. 

The  first  term  is  the  term  of  order  one  which  we  have  already  used, 

and  the  second  is  the  term  of  order  zero  which  we  have  been  seeking. 

The  point  of  3ll  of  this  work  was  to  find  the  symbol,  at  least 

to  leading  order,  of  the  operator  Z  which  represents  the  zero-order 

coupling  in  the  system  (4.15).  According  to  (4.16),  this  operator  is 
given  by 

2  =  QCQ  1  +  Q^Q  *  +  (terms  from  QHQ  *  of  order  zero  or  less)  . 

Our  results  in  (4.P)  and  (4.24)  imply  that  its  symbol  is 

(4.25)  C  +  v»MAy  +  f;(|  )  +  order  (-1)  . 

The  system  (4.15)  can  therefore  be  written  in  the  form 
3w  ? 

14.26)  -JL  =  Gwq  ♦  ZQw0  ♦  r(S£*)w  *  f'(|*)w 

+  f>(  w(l-v>))w  ♦  (order  -l)w  , 

where  2^  is  the  operator  whose  symbol  is  <p(C  +  MAy) .  We  have  split 
the  matrix  C  into  the  sum  *  (l-^)C  in  order  to  give  a  neater 
form  to  certain  formulas  which  will  appear  later. 

We  are  finally  ready  to  uncouple  the  term  of  order  zero.  To  do 
this  we  will  use  the  method  given  near  the  end  of  Section  5.2.  That 
is,  we  will  apply  an  operator  of  the  form  I  ♦  K  to  (4.26),  where  K 


has  order  -1,  and  then  make  the  corresponding  change  of  dependent 
variable.  We  will  choose  K  so  that  the  cero-order  operator  in 
the  transformed  system  has  a  diagonal  leading  symbol .  According  to 
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the  work  in  Section  3.2,  this  operator  is  KG  -  GK  +  We  there¬ 

fore  need  to  choose  K  so  that 

(4. 2"”)  -  a_a„  +  v?(C  +  MA  )  =  diagonal  matrix  . 

K  b  b  K  V 

Here  a  and  J,-  are  the  symbols  of  K  and  G,  respectively, 

N  b 

and  the  third  term  is  the  symbol  of  2^.  is  given  in  (4.15), 

A  and  C  are  given  in  (4.6),  and  M  is  given  in  (4.10). 

We  can  solve  ( 4 . 27 A  using  the  method  of  Proposition  3.2.  After 
a  certain  amount  of  labor  we  obtain 


(4.2S) 


where 


l  0 

k12 

k]  3 

U 

(1  -  f 

k21 

0 

k23 

k32 

0 

k12 

"  V? 

[(Y,!  +  Y,-)  ^ 

-cl  - 

a(ay  -  < 

:.v1] 

k13 

'T 

V  •. 

[  (' Y^l  +  Y:3)(ot  +  cl 

V’ 

k2I 

■'T 

[-a(Yj-,  +  Y-2) 

-  a 

y 

(a-c)  ] 

k_i  =  —  [:t(y12  -  y321  +  4y(a+c)] 


k23  =  4  (-Y11+V13-Y.1*Y-3Ha+c) 


k32  =  4  (Y11  +  Y13 '  V31  "  Y331 (a‘C)  ' 
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Equation  (4.2?)  does  not  impose  any  conditions  on  the  diagonal 
elemnets  of  K.  For  convenience  we  have  set  these  equal  to  zero. 

The  operator  I+K  transforms  the  system  (4.26)  into  the  form 

3w. 

(4.29)  =  GWj  +  (diagonal  operator  of  order  zerojw^  +  (order(-2))w 

O 

4m 

*  ( ■'(  +  (’{  £  <£)w  +  ^'(U)(1-^))W  , 

where  Wj  =  (I+K)Wq  =  (I+K)Qw.  The  symbol  of  K  is  given  in  (4.28), 
the  symbol  of  Q  is  given  in  (4.12),  and  the  components  of  w  are 
given  in  (4.4).  This  represents  all  of  the  uncoupling  which  we  will 
do  for  this  system. 
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4 . 2  Boundary  Conditions 

It  is  now  time  to  use  the  results  of  the  uncoupling  process 
to  derive  boundary  conditions  for  the  system  (4.1).  It  is  necessary 
to  identify  the  incoming  fast  component  for  the  partial ly  uncoupled 
system  (4.29)  and  then  find  conditions  which  suppress  this  component 
at  the  boundary  x  =  0. 

The  s>mibol  of  the  operator  G  which  appears  in  (4.29)  is  given 
in  (4.14)  and  is  equal  to 


Since  I  a  z  c  ja!.  the  second  and  third  components  of  in 

(4.29)  are  the  rapidly  moving  portions  of  the  solution.  The  second 
is  the  incoming  component,  since  a-c  <  0.  He  need  to  use  the  identity 
Wj  =  iJ  +  KlQw  to  find  an  explicit  formula  for  this  ''omponent ,  and  then 
we  need  to  use  this  formula  to  find  suitable  boundary  conditions  for 
(4.1). 

The  dependent  variable  w ^  is  given  by 


(4 . 30) 


Wj  =  (I+K)Qw 

=  "(I  +  0,1"  o  "(I  +  ~  i^M) "w  . 


Here  we  have  used  the  expression  (.4.12)  for  the  svmbcl  of  Q.  The 
quote  marks  denote  pseudo  differential  operators  having  the  stated 
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symbols,  and  the  small  circle  denotes  composition  of  operators. 

(4.301  can  be  written 

(4.311  Wj  =  "(<p  ♦  ok)"  °  "(p  +  yPM)"w  +  (  (W)  . 

In  order  to  produce  cleaner  formulas  later  on  we  have  used  the  cut¬ 
off  function  p  to  restrict  the  solution  to  the  set  7  in  the 
(u),£)  space  in  which  the  fast  waves  can  be  found. 

According  to  the  composition  law,  the  symbol  of  the  composition 
in  (4.311  is 

(p  +  "j.1  (ip  +  f  PM1  +  i  (P  +  c^l  (p  *  4  pM) 

(4.321 

+  -p  x-  (V  +  cO  -r-  (P  *  -  pM)  +  order  (-21  . 

i  7(a)  k  dy  ^ 

The  derivatives  of  p  +  c  with  respect  to  j  and  2  are  of  order  -2, 

K 

since  has  order  -1  and  we  are  ignoring  derivatives  of  p.  (4.32) 

is  therefore 

(P  +  7 , ,  1  (ip  +  4  pM)  +  order  (-21 

K  u, 

(4.331  =  p”  *  ■pc..  +  4  p“M  +  4  pc.-M  +  order  (-21 

N  L  c,  N 

=  P“  +  pc^  +  4  P”M  *  f  (  p“)  +  order  ( -21  . 

To  obt.m  the  last  line  we  used  (4.28)  to  conclude  a,.  =  (  (  i  PI  . 

i'-  ^ 

i 4. 31  I  and  '4.331  now  imply 
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'4.341  Wj  =  (£ 


^K  + 


X* 


'Ml  w  +  (  (  <p~)  *{  1 1-<£)  +  order  (-2) 


The  error  terms  in  (4.34)  can  be  ignored  in  the  system  (4.29), 

since  their  only  contributions  in  that  equation  are  error  terms 

having  the  same  order  as  terms  which  are  already  there.  In  particular, 

1 

G  (4^2)  =  U  (C)  +  rU))  0r:(-^sP2)  =  The 

3“  52  C  £2 

system  (4.29)  can  therefore  be  written 


'4.331 


=  Gz  (diagonal  of  order  zero)!  +  (order(-l)iw 


(  ~  >)w  +  (  4  >p)w  +  (w(l-v?))w  , 


wnere 


(4.36) 


*  y  v>"M)  w 


7^  is  given  in  (4.281  and  M  is  given  in  N. 10). 

The  second  component  of  z  is  the  incoming  fast  component  which 
we  need  to  suppress.  According  to  (4.361  and  the  expressions  for 
and  M,  this  component  is 


(4.3’) 


z,(x,y,t)  = 

,  f  i~  i  1 1  '  i  ■> + '  5  ^  +  ( a  -  c )  .  /  p 

*  ) i 1 T> — 


o 


ic(l  -  ^*2) 


) 


1(2) 


.  (3t+c)  (Y  .  -y,  -+T-,  -Y--1  r.y, 

,  1  ,  11  1j  el  oo  (a ) i  ,  ,  - 

(  -T~)  -  w  d—’d  z  . 

4c  . ,,  w  I 

1^(1  -  -  V?S) 

^  ) 
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Here  w ^  ' ,  w*~\  and  denote  the  components  of  w.  These  are 

(2) 

given  explicitly  in  14.4).  We  note  that  (4.37)  is  a  perturbation  of  w 
This  expression  can  be  simplified  somewhat.  The  factor 
(1  -  *  can  be  approximated  by  1  +  f(y<p).  When  this  is 

multiplied  by  (i^)  1  the  result  is  (i£)  *  +  (  (  *>) .  For  reasons 

stated  earlier,  terms  of  order  — -  can  be  omitted  without  affecting 

the  order  of  the  error  terms  in  the  system  (4.55).  We  can  therefore 
replace  (4. 3”)  with  a  new  fast  quantity 


!  '  ($)  —  -4(-~l(^v  ,+v.O  +  i  .fa-c))  i  wU) 

!  ,  j  "  1  -,  i  t  •-  O  —  } 


(4.3S) 


1“ 


Cv  - 


• (2)  1,1.,  w  , '(3)  ; 
w  -  -t  4~1  ^c)(Yll-Y13n31-Y331w  d-xlt.  . 


We  need  to  find  a  condition  which  suppresses  (4.38)  at  the  boundary 
x  =  0.  If  the  coefficients  are  independent  of  y  and  t,  then  the 
bracketed  quantity  in  the  integrand  is  the  Fourier  transform  of  (4.38), 
give  or  take  a  factor  •f",  and  we  can  accomplish  what  we  want  by 
setting  this  quantity  equal  to  zero  at  the  boundary.  If  we  do  this, 
clear  denominators,  and  then  invert  the  Fourier  transform,  we  obtain 


)t 


14.591 


~  [a(v12  +  y32)  +  av(a-c)]w(n 

c»  2 


t+c 
4  c 


1  C 


11 


v  + 
13 


31 


-  Y  -  - 1 w 


(3) 


0  for  x  =  0  . 


If  the  coefficients  depend  on  v  or  t  this  derivation  is  not 


valid.  However,  we  can  show  that  (4.39)  is  still  useful  in  this 
case.  Suppose  that  this  condition  holds,  and  write  it  in  the  simpler 
form 


(4.401 


.’WU1  j  3w(n  _  (11  „  (31 

— —  *  —  — - -  +  F,wv  +  F.w1  =  0  . 

4t  t  3v  1  a 


If  we  apply  the  operator  having  the  symbol  to  (4.401,  the  result  is 

V  '  Cl  a  _  ’•  (11  c  (31  1 

0  -  Ty  o  •  It,  W  +  (  lu)  -  +  F.  )  W  +  F  ,W 

^  a  /'=r  l  a 


(4.411 


=  "C”w' '  •  +  o  -  (  y  1  —  ♦  -ri  *  M  ~  1  +  (  |  ~  1 

(31 


(11 


•> 

V  - 

1 


"  T'-p- 
♦  <fi~  |  p-  +  ( ■(  -pr  lj  w 


ti  Ujv+itl  „ 


2  /  -j  ,  a  1  >  '  ( 1 1  '(21  3  (31- 

'r  ‘  {  (  -  1  -  ♦  —  }w'  +  W  *  t-=-  W  d-'di, 

i_  1  ,  —  i;  i 


+  (  (  j  <^“)w  +  (order(-21  lw  . 


According  to  (4.3S1  and  the  definitions  of  Fj  and  F.  implied  by 
(4.401,  the  integral  in  the  last  line  is  our  approximation  to  the  in¬ 
coming  fast  part  of  the  solution.  The  entire  last  line  is  equal  to 
zero,  so  this  fast  part  must  be  equal  to 


f  {  +  order(-21w 

at  x  =  0.  The  incoming  fast  part  is  therefore  small  compared  to  w 
for  large  frequencies  and  for  angles  near  normal  incidence. 


The  boundary  condition  (4.39)  is  written  in  terms  of  the  com¬ 


ponents  of  the  vector  w  which  appears  in  the  system  (4.3).  We 
can  use  the  definition  (4.4)  of  w  to  write  the  condition  in  terms 
of  the  variables  u,  v,  and  p  in  the  original  system  (4.1).  When 
we  do  this  the  result  is 


(4.40) 


TI  (U+P>  +  a  '  Z  [at(Y12  +  Y32)  +  Va-c)J' 


i+c 


-(ir)(Yll'Y13  +  Y3l'Y33)(u-P)=0 


for  x  =  0  . 


It  may  be  worthwhile  to  compare  (4.40)  with  the  boundary  condition 


(4.41) 


(2) 

wv  =  u  +  p  =  0 


which  we  mentioned  early  in  Section  4.1.  This  condition  was  derived 
from  the  system  (4.3),  in  which  the  coefficient  of  the  x-derivative 
is  a  diagonal  matrix.  The  newer  condition  (4.40)  is  based  on  the 
incoming  fast  variable  (4.37)  which  was  obtained  from  a  more  extensive 
uncoupling  of  the  system.  An  inspection  of  (4.57)  shows  that  this 
variable  can  be  considered  a  perturbation  of  =  u  +  p,  so  in 

some  sense  (4.40)  is  a  refinement  of  (4.41).  One  obvious  difference 
between  the  two  is  the  presence  of  derivatives  in  (4.40).  This  is  a 
result  of  the  need  to  clear  denominators  in  the  Fourier  transform 

of  the  incoming  fast  part.  The  other  difference  is  the  presence  of 

Bv  _ 

terms  which  do  not  involve  u  +  p.  The  term  a  is  a  result  of 
uncoupling  the  leading  symbol,  and  the  other  terms  in  (4.40)  are  the 


.1 
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result  of  reducing  the  coupling  caused  by  terms  of  order  zero.  The 
term  Uy(a-c)  corresponds  to  the  part  of  the  zero-order  coupling 
which  resulted  from  the  prior  uncoupling  of  the  leading  s\7nbol .  If 
we  had  not  carried  out  the  lower-order  uncoupling,  then  the  boundary 

condition  would  have  been  —  (u+p)  +  a  ^-  =  0. 

3t  r  3y 

Up  to  now  we  have  discussed  boundary  conditions  only  for  the  in¬ 
coming  fast  part  of  the  solution.  If  the  boundary  x  =  0  is  an  in¬ 
flow  boundary,  i.e.,  if  a  <  0  in  (4.31,  then  we  must  also  specify 
a  condition  for  the  slow  part.  One  possibility  is  to  use  the  system 
(4.3)  to  obtain  the  condition 

(4.42)  =  v  =  given  function,  for  x  =  0. 

Another  possibility  is  to  try  to  base  a  boundary  condition  on  the 
more  extensively  uncoupled  system  (4.35).  We  could  presumably  pre¬ 
scribe  a  value  for  the  Fourier  transform  of  the  slow  comnonent  in 
(4.55),  clear  denominators,  and  then  apply  an  inverse  transform  to 
obtain  an  inhomogeneous  boundary  condition  analogous  to  (4.40). 

The  first  approach  suggested  here  is  acceptable,  but  the  second 
one  is  not.  Our  use  of  the  cutoff  function  ^  means  that  we  have 
uncoupled  the  system  only  on  the  wedge  ?  in  the  (u),^)  space 
which  corresponds  to  rapidly  moving  waves.  This  is  clearly  no  re¬ 
striction  when  we  are  seeking  boundary  conditions  which  suppress  the 
incoming  fast  part  of  the  solution.  But  in  the  present  case  it  is 
a  major  restriction,  since  the  slow  part  of  the  solution  is  associated 
with  the  entire  (to,0  space.  The  partially  uncoupled  system  (4.55) 
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cannot  give  a  full  description  of  the  slow  part,  so  there  is  no 
point  in  trying  to  use  this  system  to  find  an  improvement  of  the 
condition  (4.42). 

We  will  therefore  prescribe  the  conditions 

-57  (u+p)  +  cx  —  -  i  [a(Yl2  +  Y32)  +  ay(a~c)lv 

(a) 

Ct+c 

(4.43)  -(-4T}  (Y11'Y13  +  Y31-Y33)(u-P)=° 

(b)  v  =  given  function,  for  x  =  0  , 

if  x  =  0  defines  an  inflow  boundary.  In  the  case  of  an  outflow 
boundary  we  will  use  only  the  first  condition. 

We  need  to  discuss  whether  these  conditions  define  a  well -posed 
initial-boundary  value  problem.  We  will  first  consider  the  special 
case  in  which  the  zero-order  terms  in  (4.43) fa)  are  not  present.  This 
would  be  the  situation  if  we  were  to  uncouple  the  leading  symbol  but 
do  nothing  about  the  zero-order  coupling  in  the  system.  In  this  case 
the  boundary  conditions  at  inflow  have  the  form 

v  =  g 

•t  ~ 

u+p  =  (u+p)t=0  -  a  (y,x)di,  for  x  =  0  , 

where  g  is  a  given  function  of  y  and  t.  Simple  energy  estimates 
for  the  system  (4.3)  show  that  this  defines  a  well-posed  problem.  The 
a  priori  estimates  for  the  solution  involve  a  time  integral  of  a 
tangential  derivative  of  the  data. 
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In  the  general  case  we  must  do  something  different  in  order  to 
properly  handle  the  incoming  fast  part  of  the  solution.  This  is 
really  no  problem,  since  in  principle  it  has  already  been  done.  In 
Section  2.6  we  described  a  process  for  estimating  the  incoming  fast 
part  for  systems  in  one  space  dimension,  and  at  the  end  of  Section 
5.2  we  indicated  that  this  process  extends  to  the  multi-dimensional 
case  with  little  modification.  These  estimates  were  derived  in  order 
to  show  that  our  boundary  conditions  are  effective  at  suppressing  the 
incoming  fast  part.  One  would  expect  that  they  also  imply  that  the 
conditions  give  well-posed  problems. 
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4 . 3  Effects  of  Orthogonal  Changes  of  Spatial  Coordinates 

The  spatial  domain  considered  in  the  previous  sections  was  the 
half-plane  for  which  x  >  0.  In  order  to  treat  slightly  more  general 
regions  we  will  now  consider  the  effects  of  linear  orthogonal  changes 
of  coordinates.  We  will  first  derive  some  general  formulas,  and  we 
will  then  use  these  formulas  to  derive  boundary  conditions  for  the 
four  sides  of  the  unit  square  0  <  x  <  1,  0  <  y  <  1.  These  conditions 

will  be  used  for  the  test  problem  which  will  be  discussed  in  the  next 
section . 

We  first  need  to  establish  some  notation.  Let  R  denote  an 
orthogonal  transformation  on  IR“,  i.e.,  either  a  rotation  or  a  flip 
of  coordinates.  This  is  illustrated  in  Figure  4.1  for  the  case  where 


R  is  a  rotation.  Suppose  that  f  is  a  function  which  is  defined 
with  respect  to  the  old  coordinates  (solid  axes),  and  let  f  be  a 
function  defined  in  the  new  coordinates  (dotted  axes'!.  From  the 
figure  we  can  see  that  if  f  is  evaluated  at  point  t  =  (x,y) ,  then 
the  numbers  plugged  into  f  must  be  given  by  ;  =  R  We  there¬ 


fore  have  the  relation 
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(4.44) 


f(z)  =  f(z)  =  f (R  z) 


for  scalar  functions.  For  vector-valued  functions  the  dependent 

variables  must  also  be  transformed.  In  this  case  the  vector 
T 

f  =  (f  ,f2)  describes  a  direction  relative  to  the  solid  axes, 
and  the  vector  f  describes  the  same  direction  relative  to  the 
dotted  axes.  The  correct  change  of  coordinates  is  therefore  given  by 


(4.45) 


f(z)  =  Rf(R  =)  - 


We  will  need  to  use  the  fact  that  tbs  'i-.dient  of  a  function  and 
the  divergence  of  a  vector  field  arf  »  in:  in  the  sense  implied 

by  (4.44)  and  (4.45).  We  will  give  %  icofs  of  these  in  order 
to  help  establish  our  notation. 

First  consider  the  gradient.  If  f(z)  =  f(R  1 z ) ,  then 


f'  (z)  =  f ’ (R_1z)R~l  , 


(  4.4b) 


(3^,3/)  =  01f,32f)R'1 


Here  the  subscripts  1  and  2  denote  differentiation  with  respect  to 
the  first  and  second  arguments,  respectively.  (4.46)  can  be  written 


(4.47) 


\32f(z)/  \  3;,f(R_1zT 


Here  we  have  used  the  fact  that  the  transformation  is  orthogonal,  i.e., 
- 1  T 

R  =  R  .  The  invariance  of  the  gradient  follows  from  the  observation 
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that  (4.4")  is  a  special  case  of  the  formula  for  changes  of  coordinates 
given  in  (4.45). 

We  now  consider  the  divergence.  If  f  and  f  are  vector-valued 
functions  such  that  f(z)  =  Rf (R  ^z),  then  f'(z)  =  Rf ' (R  *z)R  or 

Ojf.ajf)  =  ROjf.a^jR-1 . 

In  this  case  (3^f,3.,f)  is  the  Jacobian  matrix  of  f.  The  divergence 
of  f,  3jf^  +  3-,f-,,  is  equal  to  the  trace  of  this  Jacobian.  The  fact 
that  the  trace  of  a  matrix  is  invariant  under  similarity  transformations 
impl ies 

(div  f)  (z)  =  (div  f) (R_1z)  . 

We  will  now  study  the  effect  of  the  transformation  R  on  the 
system  (4.1), 


In  order  to  make  it  fit  our  notation  for  changing  coordinates,  we  will 
write  this  system  in  the  form 


(4.48)  W't  ( : ,  t ) 


‘“i 


IV. .iw  , 
n 


Y 

where  W  =  (u,v,p)  .  The  numerical  subscripts  on  u,  v,  p,  and  W 
denote  partial  derivatives.  The  3  *  2  matrix  (W^,W,)  is  the 


Jacobian  matrix  of  W  with  respect  to  the  spatial  variables.  It  can 
be  denoted  by  IV' . 

We  will  define  the  change  of  coordinates  by 


(4.491 


W(z,t) 


WCR^z.t)  . 


The  matrix  in  (4.49)  is  a 
appears  in  the  upper  left, 
necessary  to  transform  the 
ordinates  are  changed.  If 


3  *  3  matrix  in  which  the  2  *  2  matrix 
This  matrix  is  present  because  it  is 
velocity  components  when  the  spatial  co- 
(4.49)  is  inserted  into  (4.4.3),  the  result 


i  s 


In  the  first  term  on  the  right  we  used  the  chain  rule  to  evaluate  the 
Jacobian  matrix  IV' .  (4.50)  can  also  be  written 


Wt  =  (WrW2)  R 


la\  /  R  \ 

/PlU.t)\ 

' 3 1  \  if 

\p,(s.t)  * 

(4.51) 


-  c 


U1  +  V2  ' 


According  to  the  formula  (4.47)  for  the  invariance  of  the  gradient,  the 

_  -1  „  T 

second  term  on  the  right  is  c(p^(R  z,tl, p.,,0)  .  The  form  of  the 


11 


third  term  is  invariant  under  the  transformation,  since  + v, 
is  the  divergence  of  the  velocity  field.  The  system  (4.511  is 
therefore  the  same  as 


(4.521  W 


(IVj.lV.l 


(Y- -)W  , 
ij 


where 


(4 . 551 


Equation  (4.551  represents  the  transformation  of  the  velocity  field 
of  the  flow  about  which  the  system  has  been  linearized. 

Equation  (4.521  shows  that  the  form  of  the  system  (4.4S1  does 
not  change  under  orthogonal  transformations  of  spatial  coordinates. 

This  result  depends  on  the  fact  that  spatial  derivatives  appear  only 
as  gradients  and  divergences.  The  value  of  this  result  is  that  our 
previous  calculations  immediately  give  boundary  conditions  along  any 
straight  boundary.  Suppose  that  our  spatial  domain  is  the  region  to 
the  upper  right  of  the  line  l  in  Figure  4.2.  Choose  a  coordinate 
system  so  that  the  y  axis  coincides  with  Z  and  so  that  the  positive 
x  direction  points  into  the  region.  We  can  now  apply  our  earlier 
calculations  regarding  boundary  conditions  which  suppress  the  in¬ 
coming  fast  part  of  the  solution. 


1 1 


Figure  4.2 


If  this  is  an  inflow  boundary,  then  (.4.45)  implies  that  we  can  use 
37  (u+P)  +  a  ||r  -  |  [S('rl2  +  Y.,)  +  a~(5-c)]v 


(4.54) 


■(V)  (Y11-Yl3  +  Y31-Y33)(a-^  =  ° 


v  =  given  function,  for  x  =  0  . 


If  this  is  an  outflow  boundary,  then  we  would  use  only  the  first  con¬ 
dition.  The  conditions  (4.S4)  can  be  expressed  in  terms  of  the  co¬ 
efficients  and  components  of  the  original  system  (4.1)  through  the 
relations 


P  =  P 


_  1 


R 


A 
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We  will  now  use  these  general  formulas  to  derive  boundary 
conditions  for  the  sides  of  the  unit  square  0  <  x  <  1,  0  <  y  <  1. 
These  conditions  will  be  needed  for  the  test  problem  which  will  be 
discussed  in  the  next  section. 

Denote  the  sides  of  the  squares  in  the  manner  indicated  in 
Figure  4.3.  For  segment  A  we  can  use  the  conditions  (4.43)  which 


y1 

k 

D 

l 

A 

; 

i  c 

X 

00 

Figure  4.3 


were  derived  earlier  since  this  part  of  the  boundary  corresponds  to 
x  =  0.  The  inward  normal  direction  for  segment  B  is  the  positive  y 
direction,  so  for  this  segment  we  need  to  use  the  transformation  x  =  y, 
y  =  x.  The  matrix  R  is  then  given  by 


R 


B 


For  segment  C  the  transformation  is  x  =  -x,  y  =  y,  and  foi  seg¬ 
ment  D  it  is  x  =  -v,  y  =  x.  The  matrices  for  these  transforma¬ 


tions  are 
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rc=  Cl  !) and  *■>•(.!  !)■ 

respectively. 

Routine  calculations  produce  the  following  boundary  conditions: 
Segment  A  (x  =  0) : 

37  <U+P>  *  a  '  7  [a(Y12*Y32}  *  ay(a*c)]v 


(4.  55) (A) 


r— ) 

'  4c  ' 


(y11*Y 


13  y3 1 


y.3) (u-pl  =  0 


v  =  given  function. 


Segment  B  (y  =  0] : 


77  (V*P) 


3u 

3x 


t6(Y2l 


Y,J 

ol 


Sx(S-c) ]u 


(4.55} (BY 


<TT>Cv 


22  ‘  Y23  +  Y52 


Y33) (v-p)  =  0 


u  =  given  function. 


Segment  C  (x  =  1): 


-b  (-«*«?)  -  a  -  7  C-aC-Y-,  ,  *  Y,0  -  oa-ct-c)  ]v 


3t 

(4.55) (C) 


3y  c 


12  32 


(bb)  Cy„  *y,,-  y-,  -y,,)(-u-P)  «  0 


4c 


11  '15  '31  ’33 


v  =  given  function. 
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Segment  D  (y  =  1): 


) 

3t 


(-v+p) 


e  lx  '  c  t-3C-Y2i  -  Y-p 


6x(-B-c)]u 


(4.551  CL>) 


-  ( 


-k~*C 

4c 


(Y22  +  Y23 


Y32'Y33K'V-p)  =  0 


u  =  given  function. 


For  each  segment  the  first  condition  is  the  one  which  is  intended 
to  suppress  the  incoming  fast  part  of  the  solution.  The  second  con¬ 
dition  prescribes  a  value  for  the  slow  variable,  and  it  should  be 
imposed  only  when  the  segment  is  an  inflow  boundary. 


11? 


In  this  section  we  present  the  results  of  some  numerical  compu¬ 
tations  involving  the  boundary  conditions  which  have  just  been  derived. 
We  consider  the  system 


on  the  unit  square  0  <  x  <  1,  0  <  y  <  1.  This  system  is  a  special 

case  of  the  system  14.11.  We  will  use  two  different:  choices  for  the 

matrix  I  \  .  .  1 . 

ij 

We  wish  to  compare  three  types  of  boundary  conditions  for  this 
system.  The  first  of  these  is  obtained  by  diagonalicing  the  coeffi¬ 
cient  of  the  normal  derivative  and  then  defining  boundary  conditions 
in  terms  of  the  dependent  variables  in  the  new  system.  These  variables 
will  be  referred  to  as  "characteristic  variables".  This  was  discussed 
early  in  Section  4.1,  immediately  after  equation  (4.4).  For  the  four 
sides  of  the  unit  square  the  incoming  fast  characteristic  variables 
are  the  quantities  in  (4.55)  which  are  differentiated  with  respect  to 
time.  The  second  set  of  conditions  is  obtained  by  uncoupling  the 
leading  symbol  in  the  manner  described  earlier,  but  then  doing  nothing 
about  the  zero-order  coupling  in  the  system.  These  conditions  can  be 
obtained  by  deleting  the  zero-order  terms  in  the  derivative  conditions 
appearing  in  (4.55).  The  third  set  of  boundary  conditions  is  obtained 
by  also  uncoupling  the  zero-order  terms  in  the  system  of  differential 
equations.  These  are  the  conditions  (4.55). 
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We  present  two  separate  tests  of  these  conditions,  one  to  demon¬ 
strate  the  effect  of  uncoupling  the  leading  symbol ,  and  the  other  to 
demonstrate  the  effect  of  uncoupling  both  the  leading  symbol  and  the 

zero-order  term.  In  the  first  case  we  let  y.  .  =  0  for  all  i,i  and 

ij 

use  the  first  two  sets  of  boundary  conditions.  In  the  second  case  we 


use  all  three  sets  of  conditions,  and  we  let  (y .  . )  be  the  matrix 

ij 


(0  10  0  \ 

-10  0  01 

0  0  of 

In  the  computations  we  set  the  solution  equal  to  zero  when  t  =  0. 
At  the  boundary  x  =  0  we  set  v  (see  (4.55) (A))  equal  to  a  pulse 
consisting  of  half  a  sine  wave  in  t  multiplied  by  half  a  sine  wave 
in  the  tangential  variable  y.  We  use  homogeneous  conditions  on  the 
other  boundaries.  The  nonzero  part  of  the  solution  is  due  entirely 
to  the  nonzero  data  at  the  boundary  x  =  0,  so  it  is  possible  to 
study  the  influence  of  these  data  by  examining  the  size  of  the  solu¬ 
tion  in  various  parts  of  the  (x,y)  plane  at  various  times. 

In  our  computation  the  system  is  approximated  by  the  leap  frog 
difference  scheme.  The  derivative  boundary  conditions  in  (4.55)  are 
approximated  by  centered  differences  in  the  time  and  tangent  variables. 
The  outgoing  fast  characteristic  variables  are  extrapolated  at  the 
boundary  using  the  given  differential  equation.  For  this  we  use 
centered  differences  in  the  time  and  tangent  variables,  and  we  approxi¬ 
mate  the  normal  derivative  with  a  forward  difference  which  uses  a  time 


¥ 
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average  at  the  back  point.  At  an  outflow  boundary  the  slow  character¬ 
istic  variable  is  extrapolated  in  the  same  manner. 

The  boundaries  y  =  0  and  y  =  1  are  characteristic  for  the 
system  (4.56).  At  these  boundaries  we  integrate  the  slow  character¬ 
istic  variable  in  the  boundary  using  a  centered  difference  approxima¬ 
tion.  This  is  an  experiment  to  see  if  the  incoming  fast  modes  can  be 
activated  at  a  characteristic  boundary.  In  our  earlier  discussion  we 
always  assumed  that  our  boundary  was  noncharacteristic. 

The  surfaces  pictured  in  Figures  4.5  and  4.6  are  graphs  of 

•»  ~>  -i  1/  2 

(u“  +  v“  +  p')  as  functions  of  x  and  y  for  fixed  t.  The 
configuration  is  shown  in  Figure  4.4. 


Figure  4.4 

We  show  solutions  at  times  t  =  .125,  .25,  and  .375.  The  fast  mode 
entering  through  the  boundary  x  =  0  has  normal  velocity  4  since 
a  =  -1  and  c  =  -3.  Pulses  entering  on  this  mode  should  therefore 
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be  visible  near  the  nearest  boundary  (x  =  1)  in  the  graphs  for  t=  .25. 

In  Figure  4.5  we  show  the  effect  of  uncoupling  the  leading  symbol. 
In  this  case  Y—  =  0  for  all  i,j.  Figure  4.5(a)  shows  the  solution 
corresponding  to  the  boundary  conditions  defined  in  terms  of  character¬ 
istic  variables.  The  solution  in  Figure  (4.5) (b)  corresponds  to  the 
more  refined  boundary  conditions.  The  second  set  of  conditions  is 
clearly  more  effective  at  suppressing  the  incoming  fast  part  of  the 


solution. 


In  Figure  4.6  we  show  the  effect  of  uncoupling  both  the  leading 
symbol  and  the  term  of  order  zero.  In  this  case  the  matrix  (y^.) 
is  given  by  (4.57).  The  simplest  boundary  conditions  are  used  in 
part  (a).  In  part  (b)  we  use  the  boundary  conditions  obtained  by 
uncoupling  the  leading  symbol  only.  The  boundary  conditions  for  part 
(c)  are  obtained  by  uncoupling  both  the  leading  symbol  and  the  term 
of  order  zero.  The  third  set  of  conditions  is  clearly  the  most  effec¬ 


tive. 
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APPENDIX 

PROPERTIES  OF  PSEUDO  DIFFERENTIAL  OPERATORS 


In  this  appendix  we  will  define  pseudo  differential  operators 
and  state  without  proof  of  some  of  their  basic  properties.  More 
extensive  treatments  can  be  found  in  Nirenberg  [  6  ]  ,  Taylor  [  9  ]  , 
and  Treves  f  1  1  ] . 

We  must  first  establish  some  notation.  Partial  derivatives  in 
lRn  will  be  denoted  by  ?a,  where  a  =  (a, , • . . ,a  1 ,  and 


a, 

3a  =  dj1 


a. 

01  ->  1 

3  n  =  (  — ) 


JX 


1 


’(  —  ) 
^  3x  ' 
n 


The  a-  are  nonnegative  integers.  Differential  operators  can  then 
be  written  in  the  form 


P  = 


Z 


a 

a 


a 


The  coefficients  a  are  functions  on  IRn,  and  the  sum  is  taken 

a 

over  finitely  many  multi-indices  a.  We  will  allow  the  possibility 
that  P  may  act  on  vector-valued  functions.  In  that  case  the  a,, 
may  be  either  scalars  or  matrices. 

The  Fourier  transform  on  lRn  will  be  denoted  by 


u(51  = 


(2tt) 


e  lx  ^  u(x)dx  , 


]R 


where  x 


The  inverse  Fourier  transform  is  then  given  by 


ix*C  , _ 
e  u(Od£  . 

JlRn 

Differential  operators  can  be  represented  in  terms  of  the  Fourier 
transform.  For  suitable  functions  u,  we  have 

(Pu)(x)  =  Z  aa(x)3auCx) 

=  I  aa(x)3j  /  e1X‘C  G(5)dC 

=  /  elx‘C  [L  aa(x)(iC)a]  G(C)dC  . 

al  an 

Here  (i$)  1  denotes  the  product  (i£  )  •  •  U£n)  ■  The  equation 

can  be  written  in  the  form 


(A.  1) 


(Pu)Cx)  =  /  elx  '  p(x,C)G(C)dC  , 


where  p(x,£)  =  I  a^Cx)  (i£)  .  The  function  p  is  sometimes  called 
the  symbol  of  the  operator  P. 

Pseudo  differential  operators  are  obtained  by  allowing  a  larger 
class  of  symbols  to  be  used  in  (A. I).  Every  differential  operator 
is  a  pseudo  differential  operator,  but  not  vice  versa.  One  fairly 
general  symbol  class  is  the  class  S™  0  £  6  <  o  £  1,  which  was 

w  j  « 

M  OO 

introduced  by  Hormander.  This  is  defined  to  be  the  set  of  all  C 
functions  p  which  satisfy  estimates  of  the  form 


(A. 2)  3®p(x,Ol  <CK>a  6(l*  |5|)"-oM+«l8l 


x  e  k,  $  e  R 


12S 


for  all  a, 6  and  for  every  compact  subset  K  of  IRn.  The  constant 
is  allowed  to  depend  on  a.,  3,  and  K.  The  symbol  of  a  differential 
operator  of  order  m  having  smooth  coefficients  clearly  belongs  to 
the  class  S^1  q.  This  class  will  also  be  denoted  by  Sm.  For  the 
operators  considered  in  this  paper  we  always  have  p  =  1  and  5=0. 
In  general,  the  number  m  appearing  in  (A. 2)  is  called  the  order  of 
the  operator  P  whose  symbol  is  p.  The  order  need  not  be  positive, 
and  it  need  not  be  an  integer. 

If  u  e  C*,  then  Pu  e  c”.  It  is  possible  to  extend  P  so 
that  Pu  is  defined  for  any  distribution  u  having  compact  support. 
In  this  case  Pu  is  a  distribution. 

A  useful  concept  is  that  of  an  asymptotic  expansion  of  a  symbol. 

Suppose  that  {m_. ^  is  a  sequence  of  real  numbers  such  that 

m.  >  m.  .  for  all  i  and  m.  -+■  -<*  as  i  -*■  ».  Let  ;p.’1  be  a 
ii+l  •  i  •  i 

HI; 

sequence  of  symbols  such  that  p.  e  S  -  for  each  i.  A  symbol  is 
said  to  be  an  asymptotic  sum  of  the  p.. ,  written 


00 


k  mVr  + 1 

provided  p  -  E.  _  p.  e  S  for  dll  k.  That  is,  the  error  in 

j  =  0 

each  partial  sum  must  have  the  same  order  as  the  first  term  omitted 
from  the  partial  sum.  This  concept  is  analogous  to  the  usual  concept 
of  asymptotic  expansion.  In  fact,  if  a  function  p(51  of  one 
variable  has  an  asymptotic  expansion 
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“  a . 

P(Q  ~  l  4  as  c  +  “ 

j=o 

in  the  usual  sense,  then  this  expansion  is  also  asymptotic  in  the 
sense  described  above. 

Pseudo  differential  operators  can  be  composed.  Let  P  and  Q 
be  operators  with  symbols  p(x,£)  and  q(x,£)  ,  respectively,  and 
suppose  that  q  has  compact  support  in  x.  The  composition  P(Qu) 
is  then  well-defined  and  is  given  by  a  pseudo  differential  operator 
whose  symbol  has  the  asymptotic  expansion 

(A.  3)  ap  -  l  —.4 -  3rP(x,-)3*q(x,‘)  . 

lal>0  i  a  a!  ^  x 

The  sum  is  taken  over  all  multi-indices  a  =  (:<j . a  )  having  non¬ 

negative  components.  The  order  of  a  is  given  by  Lx!  =  la.,  and 
the  factorial  a!  denotes  the  product  •  ...  *0^!. 

It  follows  from  (A. 2)  that  when  a  symbol  is  differentiated  with 
respect  to  the  result  is  a  symbol  of  lower  order.  This  implies 

that  the  leading  order  term  in  (A. 3)  corresponds  to  j.x[  =  0  3nd  is 
equal  to  p(x, £)q (x, •  The  symbol  of  the  product  of  two  operators 
is  therefore  equal  to  the  product  of  their  symbols,  up  to  certain 
terms  of  lower  order. 

This  makes  sense  when  we  consider  the  special  case  of  differential 

ot  3 

operators.  The  composition  of  two  operators  a(x)3  and  b(x)3 
is  equal  to 
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^  ^  p 

a3  (b3  )  =  ab3  +  lower  order  terms  involving  derivatives  of  b. 


In  this  case  the  composition  law  can  be  derived  using  Leibniz'  rule. 
For  general  pseudo  differential  operators  the  derivation  is  much  more 
complicated. 

It  is  sometimes  necessary  to  discuss  adjoints  of  pseudo  differen¬ 
tial  operators.  The  adjoint  of  an  operator  P  is  a  pseudo  differen¬ 
tial  operator  P*  whose  symbol  has  the  asymptotic  expansion 


l 

al  >  0 


353ip*(x’c)  • 


The  leading  order  term  in  this  expansion  is  the  adjoint  of  the  symbol 
of  P. 
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